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Abstract 


Verification  and  validation  methodology  is  presented  for  CFD  simulation  results  from 
an  already  developed  RANS  CFD  code  applied  for  specified  objectives,  geometry, 
conditions,  and  available  benchmark  information.  Concepts  and  definitions  are  provided 
for  errors  and  uncertainties  and  verification  and  validation.  The  simulation  error  and 
uncertainty  equations  are  derived  with  modeling  and  numerical  errors  being  additive  and 
modeling  and  numerical  uncertainties  combining  by  root-sum-square.  The  concepts  and 
definitions  provide  the  mathematical  framework  for  the  verification  and  validation 
methodology. 

Verification  is  defined  as  a  process  for  assessing  numerical  uncertainty  and,  when 
conditions  permit,  estimating  the  sign  and  magnitude  of  the  numerical  error  itself  and  the 
uncertainty  in  that  error  estimate.  Iterative  and  parameter  convergence  studies  are 
conducted  using  multiple  solutions  with  systematic  parameter  refinement  to  estimate 
numerical  errors  and  uncertainties.  Three  convergence  conditions  are  possible:  (i) 
monotonic  convergence;  (ii)  oscillatory  convergence;  and  (iii)  divergence.  For  condition 
(i),  generalized  Richardson  extrapolation  for  J  input  parameters  and  use  of  correction 
factors  to  account  for  the  effects  of  higher-order  terms  and  defining  and  estimating  errors 
and  uncertainties  is  used.  For  condition  (ii),  the  upper  and  lower  bounds  of  the  solution 
oscillation  are  used  to  estimate  uncertainties.  For  condition  (iii),  errors  and  uncertainties 
can  not  be  estimated. 

Validation  is  defined  as  a  process  for  assessing  modeling  uncertainty  by  using 
benchmark  experimental  data  and,  when  conditions  permit,  estimating  the  sign  and 
magnitude  of  the  modeling  error  itself.  The  comparison  error  (difference  between  data  and 
simulation  values)  and  validation  uncertainty  (combination  of  uncertainties  in  data  and 
portion  of  simulation  uncertainties  that  can  be  estimated)  are  used  in  this  process. 

An  example  is  provided  for  a  RANS  CFD  code  and  results  for  steady  flow  for  a 
cargo/container  ship. 
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1.  Introduction 

Discussion  and  methodology  for  estimating  errors  and  uncertainties  in 
computational  fluid  dynamics  (CFD)  simulations  has  reached  a  certain  level  of  maturity 
with  increased  attention  and  recent  progress  on  common  concepts  and  terminology 
(AIAA,  1998),  advocacy  and  detailed  methodology  (Roache,  1998),  and  numerous  case 
studies  (e.g.,  Mehta,  1998).  Progress  has  been  accelerated  in  response  to  the  urgent  need 
for  achieving  consensus  on  concepts  and  terminology  and  useful  methodology,  as  CFD  is 
applied  to  increasingly  complex  geometry  and  physics  and  integrated  into  the  engineering 
design  process.  Such  consensus  is  required  to  realize  the  goals  of  simulation-based  design 
and  other  uses  of  CFD  such  as  simulating  flows  for  which  experiments  are  difficult  (e.g., 
full-scale  Reynolds  numbers,  hypersonic  flows,  off-design  conditions).  In  spite  of  the 
progress  and  urgency,  the  various  viewpoints  have  not  converged  and  current 
methodology  falls  short  of  providing  practical  procedures  and  methodology  for  estimating 
errors  and  uncertainties  in  CFD  simulations. 

The  present  work  provides  a  pragmatic  approach  for  estimating  errors  and 
uncertainties  in  CFD  simulations.  Previous  work  on  verification  (Stem  et  al.,  1996)  is 
extended  and  put  on  a  more  rigorous  foundation  and  combined  with  subsequent  work  on 
validation  (Coleman  and  Stem,  1997)  [hereafter  referred  to  as  C&S]  thereby  providing  the 
framework  for  overall  procedures  and  methodology.  The  philosophy  is  strongly 
influenced  by  experimental  fluid  dynamics  (EFD)  uncertainty  analysis  (Coleman  and 
Steele,  1999),  which  has  been  standardized.  Hopefully,  CFD  verification  and  validation 
procedures  and  methodology  can  reach  a  similar  level  of  maturity  and  user  variability  can 
reach  similar  low  levels,  as  for  EFD. 

The  work  is  part  of  a  larger  program  (Rood,  1996)  for  developing  and 
implementing  a  strategy  for  verification  and  validation  of  Reynolds-averaged  Navier- 
Stokes  (RANS)  ship  hydrodynamics  CFD  codes.  The  program  includes  complementary 
CFD  and  EFD  towing-tank  investigations  and  considers  errors  and  uncertainties  in  both 
the  simulations  and  the  data  in  assessing  the  success  of  the  verification  and  validation 
efforts.  The  work  also  benefited  from  collaboration  with  the  21st  and  22nd  International 
Towing  Tank  Resistance  Committees  (ITTC,  1996  and  1999). 

The  focus  is  on  verification  and  validation  procedures  and  methodology  for  CFD 
simulation  results  from  an  already  developed  CFD  code  applied  for  specified  objectives, 
geometry,  conditions,  and  available  benchmark  information.  The  procedures  and 
methodology  were  developed  considering  RANS  CFD  codes,  but  should  be  applicable  to 
a  fairly  broad  range  of  codes  such  as  boundary-element  methods  and  certain  aspects  of 
large-eddy  and  direct  numerical  simulations. 

The  present  work  differs  in  many  respects  from  recent  literature.  The  presentation 
is  relatively  succinct  with  intention  for  use  for  practical  applications  (i.e.,  industrial  CFD) 
for  which  numerical  errors  and  uncertainties  can  not  be  considered  negligible  or 
overlooked.  The  definitions  of  errors  and  uncertainties  and  verification  and  validation  that 
are  used  in  any  approach  need  to  be  clearly  stated.  Table  1  summarizes  the  present 
definitions  along  with  those  given  by  the  AIAA  (1998)  and  Roache  (1998)  for 
comparison.  The  present  and  Roache  (1998)  definitions  for  errors  and  uncertainties  are 
consistent  with  those  used  for  EFD.  The  AIAA  (1998)  definitions  are  from  an  information 
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theory  perspective  and  differ  from  those  used  in  EFD,  but  are  not  contradictory  to  the 
present  definitions.  The  present  definitions  for  verification  and  validation  are  closely  tied 
to  the  present  definitions  of  errors  and  uncertainties  and  equations  derived  for  simulation 
errors  and  uncertainties.  The  Roache  (1998)  and  AIAA  (1998)  definitions  are  broader, 
but  not  contradictory  to  the  present  definitions.  The  present  approach  includes  both  the 
situations  (1)  of  estimating  errors  and  the  uncertainty  of  those  estimates  and  (2)  of 
estimating  uncertainties  only.  Richardson  extrapolation  (RE)  is  used  for  verification, 
which  is  not  new;  however,  the  present  generalizations  for  J  input  parameters  and  use  of 
correction  factors  to  account  for  the  effects  of  higher-order  terms  and  in  defining  and 
estimating  errors  and  uncertainties  constitute  a  new  approach.  The  use  of  quantitative 
estimates  for  errors  and  the  use  of  uncertainties  for  those  estimates  also  constitute  a  new 
approach  in  verification  and  validation. 


2,  Verification  and  Validation  Procedures 

The  overall  CFD  verification  and  validation  procedures  can  be  conveniently  grouped  in 
four  consecutive  steps:  (1)  preparation;  (2)  verification;  (3)  validation;  and  (4) 
documentation. 

Preparation.  The  1st  step  is  preparation,  which  involves  selection  of  the  CFD  code 
and  specification  of  objectives,  geometry,  conditions,  and  available  benchmark 
information.  The  objectives  might  be  prediction  of  certain  variables  at  certain  levels  of 
validation  (e.g.,  programmatic  validation  requirements  U  .).  The  variables  can  either  be 

integral  (e.g.,  resistance)  or  point  (e.g.,  mean  velocities  and  turbulent  Reynolds  stresses) 
values  and  the  programmatic  validation  requirements  may  be  different  for  each  variable. 

Verification.  The  2nd  step  is  verification,  which  is  defined  as  a  process  for  assessing 
simulation  numerical  uncertainty  USN  and,  when  conditions  permit,  estimating  the  sign 

and  magnitude  5  stv  of  the  simulation  numerical  error  itself  and  the  uncertainty  in  that 
error  estimate  (referred  to  as  the  corrected  simulation  numerical  uncertainty  UScN). 

Iterative  and  input  parameter  convergence  studies  are  conducted  using  multiple  solutions 
with  systematic  parameter,  as  described  in  Section  3.2. 

Validation.  The  3rd  step  is  validation,  which  is  defined  as  a  process  for  assessing 
simulation  modeling  uncertainty  U SM  by  using  benchmark  experimental  data  and,  when 

conditions  permit,  estimating  the  sign  and  magnitude  of  the  simulation  modeling  error  SSM 
itself.  The  comparison  error  E  (difference  between  data  D  and  simulation  S  values)  and 
validation  uncertainty  Uv  (combination  of  uncertainties  in  data  and  portion  of  simulation 
uncertainties  that  can  be  estimated)  are  used,  as  described  in  Section  3.3. 

Documentation.  The  4th  step  is  documentation,  which  is  detailed  presentation  of  the 
CFD  code  (equations,  initial  and  boundary  conditions,  modeling,  and  numerical  methods), 
objectives,  geometry,  conditions,  verification,  validation,  and  analysis. 
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3.  Verification  and  Validation  Methodology 

Verification  (Section  3.2)  and  validation  (Section  3.3)  methodology  is  presented  for 
CFD  simulation  results  from  an  already  developed  CFD  code  applied  for  specified 
objectives,  geometry,  conditions,  and  available  benchmark  information.  Section  3.1 
discusses  concepts  and  definitions  for  errors  and  uncertainties  and  verification  and 
validation,  which  provide  the  mathematical  framework  for  the  verification  and  validation 
methodology.  Analytical  benchmarks  can  be  defined  as  the  truth  and  are  useful  in 
development  and  confirmation  of  verification  procedures  and  methodology  and  in  code 
development,  but  can  not  be  used  for  validation  and  are  restricted  to  simple  equations. 
Results  from  the  use  of  analytical  benchmarks  are  provided  in  Appendix  C. 


3,1  Concepts  and  Definitions 

Accuracy  indicates  the  closeness  of  agreement  between  a  simulation/experimental 
value  of  a  quantity  and  its  true  value.  Error  8  is  the  difference  between  a  simulation  value 
or  an  experimental  value  and  the  truth.  Accuracy  increases  as  error  approaches  zero.  The 
true  values  of  simulation/experimental  quantities  are  rarely  known.  Thus,  errors  must  be 
estimated.  An  uncertainty  U  is  an  estimate  of  an  error  such  that  the  interval  ±  U  contains 
the  true  value  of  8  95  times  out  of  100.  An  uncertainty  interval  thus  indicates  the  range  of 
likely  magnitudes  of  8  but  no  information  about  its  sign. 

For  simulations,  under  certain  conditions,  errors  can  be  estimated  including  both  sign 
and  magnitude  (referred  to  as  an  error  estimate  8*).  Then,  the  uncertainty  considered  is 
that  corresponding  to  the  error  in  8*.  When  8  is  estimated,  it  can  be  used  to  obtain  a 
corrected  value  of  the  variable  of  interest. 

Sources  of  errors  and  uncertainties  in  results  from  simulations  can  be  divided  into  two 
distinct  sources:  modeling  and  numerical.  Modeling  errors  and  uncertainties  are  due  to 
assumptions  and  approximations  in  the  mathematical  representation  of  the  physical 
problem  (such  as  geometry,  mathematical  equation,  coordinate  transformation,  boundary 
conditions,  turbulence  models,  etc.)  and  incorporation  of  previous  data  (such  as  fluid 
properties)  into  the  model.  Numerical  errors  and  uncertainties  are  due  to  numerical 
solution  of  the  mathematical  equations  (such  as  discretization,  artificial  dissipation, 
incomplete  iterative  and  grid  convergence,  lack  of  conservation  of  mass,  momentum,  and 
energy,  internal  and  external  boundary  non-continuity,  computer  round-off,  etc.).  The 
present  work  assumes  that  all  correlations  among  errors  are  zero,  which  is  doubtless  not 
true  in  all  cases,  but  the  effects  are  assumed  negligible  for  the  present  analyses. 

The  simulation  error  8  s  is  defined  as  the  difference  between  a  simulation  result  S  and 
the  truth  T.  In  considering  the  development  and  execution  of  a  CFD  code,  it  can  be 
postulated  that  8  s  is  comprised  of  the  addition  of  modeling  and  numerical  errors 

s5=s-r=sw+s„  (i) 

A  derivation  of  the  simulation  error  equation  (1)  is  provided  in  Appendix  A.  The 
uncertainty  equation  corresponding  to  error  equation  (1)  is 

U2  =U2  +U2  (2) 
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where  Us  is  the  uncertainty  in  the  simulation  and  USM  and  USN  are  the  simulation 
modeling  and  numerical  uncertainties. 

For  certain  conditions,  the  numerical  error  8  SN  can  be  considered  as 


'  SN 


’ SN 


+  £ 


SN 


(3) 


where  8*SN  is  an  estimate  of  the  sign  and  magnitude  of  8  s  v  and  £Sn  is  the  error  in  that 

estimate  (and  is  estimated  as  an  uncertainty  since  only  a  range  bounding  its  magnitude  and 
not  its  sign  can  be  estimated).  The  corrected  simulation  value  S(  is  defined  by 

sc=s-s;,  (4) 


with  error  equation 

The  uncertainty  equation  corresponding  to  error  equation  (5)  is 


^sr~^c  T  —  §sm+£sn 


jj 2  =U2  +U2 

u  Sr  u  SM  SrN 


(5) 

(6) 


where  U  Sc  is  the  uncertainty  in  the  corrected  simulation  and  U  s  N  is  the  uncertainty 
estimate  for  zSN. 

Debate  on  verification  and  validation  has  included  discussion  on  whether  errors  such 
as  8  SN  are  deterministic  vs.  stochastic  and  thus  how  they  should  be  treated  in  uncertainty 
analysis  was  unclear.  In  the  approach  given  by  equations.  (3)-(6),  a  deterministic  estimate 
8  *N  of  8  SN  and  consideration  of  the  error  £Sn  in  that  estimate  are  used.  The  approach  is 

analogous  to  that  in  EFD  when  an  asymmetric  systematic  uncertainty  is  “zero-centered” 
by  inclusion  of  a  model  for  the  systematic  error  in  the  data  reduction  equation  and  then  the 
uncertainty  considered  is  that  associated  with  the  model  (Coleman  and  Steele,  1999).  In 
the  “uncorrected”  approach  given  by  equations  (l)-(2),  any  particular  8  SN  is  considered  as 

a  single  realization  from  some  parent  population  of  8  SN  's  and  the  uncertainty  USN  is 
interpreted  accordingly  in  analogy  to  the  estimation  of  uncertainties  in  EFD  (with  a  similar 
argument  for  £Sn  and  Us  N). 


3.2  Verification 

For  many  CFD  codes,  the  most  important  numerical  errors  and  uncertainties  are  due 
to  use  of  iterative  solution  methods  and  specification  of  various  input  parameters  such  as 
spatial  and  time  step  sizes  and  other  parameters  (e.g.,  artificial  dissipation).  The  errors 
and  uncertainties  are  highly  dependent  on  the  specific  application  (geometry  and 
conditions). 

The  errors  due  to  specification  of  input  parameters  are  decomposed  into  error 
contributions  from  iteration  number  8  7 ,  grid  size  8  G ,  time  step  8  T ,  and  other  parameters 
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8P,  which  gives  the  following  expressions  for  the  simulation  numerical  error  and 
uncertainty 

&sn  =5/+5G+5r+5p  =5/+^8i  (7) 

M 

U2SN=Uf  +U2g+U2t  +U2p  =U2  +'£u2  (8) 

M 

Similarly,  error  estimates  8*  can  be  decomposed  as 

which  gives  the  following  expressions  for  the  corrected  simulation  and  corrected 
simulation  numerical  uncertainty 

5c =s-(8;+£s;)=r+sa,  +e„  do) 


y=i 


Verification  is  based  on  equation  (10),  which  is  put  in  the  form 

s  =  sc+(  s;+£s;)  (i2) 

M 

Equation  (12)  expresses  S  as  the  corrected  simulation  value  Sc  plus  numerical  errors.  Sc  is 
also  referred  to  as  a  numerical  benchmark  since  it  is  equal,  as  shown  by  equation  (10),  to 
the  truth  plus  simulation  modeling  error  and  presumable  small  error  £Sn  in  the  estimate  of 
the  numerical  error  5  *v .  Power-series  expansions  for  each  input  parameter  and  multiple 

solutions  are  used  to  obtain  estimates  for  the  8  * ’s  in  equation  (12).  For  this  approach  to 
be  useful,  8  *  must  be  accurately  estimated  or  be  negligible  for  each  solution. 


3.2.1  Convergence  Studies 

Iterative  and  parameter  convergence  studies  are  conducted  using  multiple  (m) 
solutions  and  systematic  parameter  refinement  by  varying  the  kth  input  parameter 
Axk  while  holding  all  other  parameters  constant.  The  present  work  assumes  input 

parameters  can  be  expressed  such  that  the  finest  resolution  corresponds  to  the  limit  of 
infinitely  small  parameter  values.  Many  common  input  parameters  are  of  this  form,  e.g., 
grid  spacing,  time  step,  and  artificial  dissipation.  Additionally,  a  uniform  parameter 
refinement  ratio  rk  =  AxkjAxki  =AxkjAxk^  =  Axkm/Axkm between  solutions  is  assumed. 
The  use  of  uniform  parameter  refinement  ratio  is  not  required  (Roache,  1998);  however,  it 
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simplifies  the  analysis  and  in  the  authors  experience  use  of  non-uniform  parameter 
refinement  ratio  is  not  needed. 

Careful  consideration  should  be  given  to  selection  of  uniform  parameter  refinement 
ratio.  The  most  appropriate  values  for  industrial  CFD  are  not  yet  fully  established.  Small 
values  (i.e.,  very  close  to  one)  are  undesirable  since  solution  changes  will  be  small  and 
sensitivity  to  input  parameter  may  be  difficult  to  identify  compared  to  iterative  errors. 
Large  values  alleviate  this  problem;  however,  they  also  may  be  undesirable  since  the  finest 
step  size  may  be  prohibitively  large  if  the  coarsest  step  size  is  designed  for  sufficient 
resolution  such  that  similar  physics  are  resolved  for  all  m  solutions.  Also,  similarly  as  for 
small  values,  solution  changes  for  the  finest  step  size  may  be  difficult  to  identify  compared 
to  iterative  errors  since  iterative  convergence  is  more  difficult  for  small  step  size.  Another 
issue  is  that  for  parameter  refinement  ratio  other  than  rk  =  2,  interpolation  to  a  common 

location  is  required  to  compute  solution  changes,  which  introduces  interpolation  errors. 
Roache  (1998)  discusses  methods  for  evaluating  interpolation  errors.  However,  for 

industrial  CFD,  rk  =  2  may  often  be  too  large.  A  good  alternative  may  be  rk  =  y[2  ,  as  it 

provides  fairly  large  parameter  refinement  ratio  and  at  least  enables  prolongation  of  the 
coarse-parameter  solution  as  an  initial  guess  for  the  fine-parameter  solution. 

Equation  (12)  is  written  for  the  kth  parameter  and  mth  solution  as 

s„.  =sc+K.  +K+  XX.  <13> 

Iterative  convergence  must  be  assessed  and  Sk  corrected  for  iterative  errors  prior  to 

evaluation  of  parameter  convergence  since  the  level  of  iterative  convergence  may  not  be 
the  same  for  all  m  solutions  used  in  the  parameter  convergence  studies.  Methods  for 
estimating  U  k  or  8  *  and  U l(  are  described  in  Section  3.2.2.  With  8*  evaluated,  Sk  is 

corrected  for  iterative  errors  as 

s,„=s>„- k  »«c+*;+  xx  <w) 

M -j*k 

Equation  (13)  shows  that  iterative  errors  8  *  must  be  accurately  estimated  or  negligible  in 

comparison  to  8  *k  for  accurate  convergence  studies  and  that  they  should  be  considered 
within  the  context  of  convergence  studies  for  each  input  parameter. 

Sk  can  be  calculated  for  both  integral  (e.g.,  resistance  coefficients)  and  point  (e.g., 

surface  pressure,  wall-shear  stress,  and  velocity)  variables.  Sk  can  be  presented  as  an 

absolute  quantity  (i.e.,  non-normalized)  or  normalized  with  the  solution  as  a  percentage 
change;  however,  if  the  solution  value  is  small,  a  more  appropriate  normalization  may  be 
the  range  of  the  solution. 

Convergence  studies  require  a  minimum  of  m=3  solutions  to  evaluate  convergence 
with  respect  to  input  parameter.  Note  that  m=2  is  inadequate,  as  it  only  indicates 
sensitivity  and  not  convergence,  and  that  m>3  may  be  required.  Consider  the  situation  for 
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3  solutions  corresponding  to  fine  Sk  ,  medium  St  ,  and  coarse  Sk  values  for  the  kth  input 

parameter.  Solution  changes  £  for  medium-fine  and  coarse-medium  solutions  and  their 
ratio  Rk  are  defined  by 

82it  =  ft*  - 

e32  k=Sk]~Sk2  (15) 

Rk  =£21t  /e32t 

Three  convergence  conditions  are  possible: 

(i)  Monotonic  convergence:  0  <  Rk  <1 

(ii)  Oscillatory  convergence:  Rk  <  01  (16) 

(iii)  Divergence:  Rk  >  1 

For  monotonic  convergence  (i),  generalized  RE  is  used  to  estimate  Uk  or  6  *  and 
Uk  .  Methods  for  estimating  errors  and  uncertainties  for  condition  (i)  are  described  in 
Section  3.2.3. 

For  oscillatory  convergence  (ii),  the  solutions  exhibit  oscillations,  which  may  be 
erroneously  identified  as  condition  (i)  or  (iii).  This  is  apparent  if  one  considers  evaluating 
convergence  condition  from  three  points  on  a  sinusoidal  curve  (Coleman  et  al.,  1999). 
Depending  on  where  the  three  points  fall  on  the  curve,  the  condition  could  be  incorrectly 
diagnosed  as  either  monotonic  convergence  or  divergence.  Methods  discussed  here  for 
estimating  uncertainties  Uk  for  condition  (ii)  require  more  than  m= 3  solutions  and  are 
described  in  Section  3.2.4. 

For  divergence  (iii),  the  solutions  diverge  and  errors  and  uncertainties  can  not  be 
estimated.  Additional  remarks  are  given  in  Section  3.2.5. 

Determination  of  the  convergence  ratio  Rk  for  point  variables  can  be  problematic  since 
solution  changes  e21  and  e32t  can  both  go  to  zero  (e.g.,  in  regions  where  the  solution 

contains  an  inflection  point).  In  this  case,  the  ratio  becomes  ill  conditioned.  However,  the 
convergence  ratio  can  be  used  in  regions  where  the  solution  changes  are  both  non-zero 
(e.g.,  local  solution  maximums  or  minimums).  Another  approach  is  to  use  a  global 
convergence  ratio  Rk,  which  overcomes  ill  conditioning,  based  on  the  F2  norm  of  the 

solution  changes,  i.e.,  (Rk )  =  e2ik  J  ^  yik  2  ■  <  >  is  used to  denote  an  averaged  value  and 
rv  iI/2 

||e||^  =  ^  £  2  denotes  the  F2  norm  of  solution  change  over  the  N  points  in  the  region 
L  i=i 

of  interest.  Caution  should  be  exercised  when  defining  the  convergence  ratio  from  the 
ratio  of  the  F2  norm  of  solution  changes  because  the  oscillatory  condition  (Rk  <  1)  cannot 

1  As  discussed  in  the  text  that  follows,  0  <  Rk  <  1  and  Rk  >  1  may  also  occur  for  the  oscillatory 
condition. 
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be  diagnosed  since  (i?; )  will  always  be  greater  than  zero.  Local  values  of  Rk  at  solution 

maximums  or  minimums  should  also  be  examined  to  confirm  the  convergence  condition 
based  on  an  L2  norm  definition. 


3.2.2  Iterative  Convergence 

Iterative  convergence  must  be  assessed  and  simulation  results  Sk  corrected  for 

iterative  errors  prior  to  evaluation  of  parameter  convergence  since  the  level  of  iterative 
convergence  may  not  be  the  same  for  all  m  solutions  used  in  the  parameter  convergence 
studies.  Methods  for  estimating  U  k  or  8*  and  Uj  are  described  in  this  section.  The 

methods  are  applicable  to  both  integral  and  point  variables.  For  point  variables,  an  L2 
norm  over  all  grid  points  is  often  used  as  a  global  metric.  There  are  many  integral  and 
point  variables  that  can  be  monitored  to  establish  iterative  stopping  criteria;  however, 
present  discussion  is  specifically  within  the  context  of  evaluating  U ,  or  8  *  and  U ,  for 

use  in  the  parameter  convergence  study  for  Sk  .  Further  work  is  needed  on  assessing 

iterative  errors  and  their  role  in  parameter  convergence  studies  and  for  assessing  iterative 
errors  and  uncertainties  for  unsteady  flows. 

Typical  CFD  solution  techniques  for  obtaining  steady  state  solutions  involve  beginning 
with  an  initial  guess  and  performing  time  marching  or  iteration  until  a  steady  state  solution 
is  achieved.  For  time-accurate  calculations  using  implicit  methods,  convergence  of  the 
solution  is  required  at  each  time  step.  Care  must  be  exercised  in  evaluating  iterative 
convergence  based  solely  on  solution  residuals,  i.e.,  change  in  solution  from  iteration  to 
iteration.  Small  time  steps  and/or  relaxation  parameters  can  result  in  small  solution 
residuals  while  iterative  error  can  be  large  (Ferziger  and  Peric,  1997).  If  Sk  is  a  primary 

dependent  variable,  an  alternative  approach  that  removes  this  problem  is  to  use  the 
residual  imbalance  of  the  discretized  equations  (i.e.,  the  difference  in  the  left-  and  right- 
hand  sides)  as  a  measure  of  convergence;  since,  the  iterative  error  satisfies  the  same 
equation  as  this  residual  imbalance. 

The  number  of  order  magnitude  drop  and  final  level  of  solution  residual  (or  residual 
imbalance)  can  be  used  to  determine  stopping  criteria  for  iterative  solution  techniques. 
Iterative  convergence  to  machine  zero  is  desirable,  but  for  complex  geometry  and 
conditions  it  is  often  not  possible.  Three  or  four  orders  of  magnitude  drop  in  solution 
residual  to  a  level  of  10"4  is  more  likely  for  these  cases.  Methods  for  estimation  of  iterative 
errors  and  uncertainties  can  be  based  on  graphical,  as  discussed  below,  or  theoretical 
approaches  and  are  dependent  on  the  type  of  iterative  convergence:  (a)  oscillatory;  (b) 
convergent;  or  (c)  mixed  oscillatory /convergent. 

For  oscillatory  iterative  convergence  (a),  the  deviation  of  the  variable  from  its  mean 
value  provides  estimates  of  the  iterative  uncertainty  based  on  the  range  of  the  maximum 
Sv  and  minimum  SL  values 


\(SV-SL) 


(17) 
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For  convergent  iterative  convergence  (b),  a  curve-fit  of  an  exponential  function 
can  be  used  to  estimate  U ,  or  8  *  and  U 7  as  the  difference  between  the  value  and  the 

exponential  function  from  a  curve  fit  for  large  iteration  number  CF x 

Uj  =  IS-CFJ 

(18) 

s;  =s-cf„,u,  =0 

1k„,  lC 


For  mixed  convergent/oscillatory  iterative  convergence  (c),  the  amplitude  of  the 
solution  envelope  decreases  as  the  iteration  number  increases,  the  solution  envelope  is 
used  to  define  the  maximum  Sv  and  minimum  SL  values  in  the  Ith  iteration,  and  to 

estimate  U l  or  8  *  and  U , 


1 

U,= 

2  ^V-SL) 

s;  =s~- 

A,  2 

■C sv-sL),u 

(19) 


An  increase  in  the  amplitude  of  the  solution  envelope  as  the  iteration  number 
increases  indicates  that  the  solution  is  divergent. 

Estimates  of  the  iterative  error  based  on  theoretical  approaches  are  presented  in 
Ferziger  and  Peric  (1997)  and  involve  estimation  of  the  principal  eigenvalue  of  the 
iteration  matrix.  The  approach  is  relatively  straightforward  when  the  eigenvalue  is  real  and 
the  solution  is  convergent.  For  cases  in  which  the  principal  eigenvalue  is  complex  and  the 
solution  is  oscillatory  or  mixed,  the  estimation  is  not  as  straightforward  and  additional 
assumptions  are  required. 


3.2.3  Monotonic  Convergence:  Generalized  Richardson  Extrapolation 

For  monotonic  convergence,  i.e.,  condition  (i)  in  equation  (16),  generalized  RE  is 
used  to  estimate  Uk  or  8  *  and  Uk  .  RE  is  generalized  for  J  input  parameters  and  use  of 

correction  factors  to  account  for  the  effects  of  higher-order  terms  and  defining  and 
estimating  errors  and  uncertainties,  as  summarized  in  the  following.  Appendix  B  provides 
a  detailed  description. 

Generalized  RE  begins  with  equation  (14).  The  error  terms  on  the  right-hand-side 
of  equation  (14)  are  of  known  form  (i.e.,  power  series  expansion  in  Axk )  based  on 

analysis  of  the  modified  (A. 6)  and  numerical  error  (A. 9)  equations,  as  shown  in  Appendix 
A  equation  (A.  12),  which  is  written  below  as  a  finite  sum  (i.e.,  error  estimate)  and  for  the 
kth  parameter  and  mth  solution 

K.  =X>V„>'’rrf’  (20) 

1=1 

n  =  number  of  terms  retained  in  the  power  series,  powers  /.(''correspond  to  order  of 
accuracy  (for  the  ith  term),  and  g(k  are  referred  to  as  “grid”  functions  which  are  a 
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function  of  various  orders  and  combinations  of  derivatives  of  S  with  respect  to  xk. 
Substituting  equation  (20)  into  equation  (14)  results  in 

4  =sc+£(A*t,  Is;,  (21) 

i= 1 

Subtraction  of  multiple  solutions  where  input  parameter  Any.  is  uniformly  refined 
eliminates  the  SJ  terms  in  equation  (21)  since  SJ  is  independent  of  Ax,,  and  provides 

equations  for  Sc,  pk] ,  and  gk] .  This  assumes  p(k  and  g‘tn  are  also  independent  of  Ary  . 

Since  each  term  (i)  contains  2  unknowns,  m=2n+l  solutions  are  required  to  estimate  the 
numerical  benchmark  Sc  and  the  first  n  terms  in  the  expansion  in  equation  (21)  (i.e.,  for 
n=l,  777=3  and  for  n= 2,  m= 5,  etc).  The  accuracy  of  the  estimates  depends  on  how  many 
terms  are  retained  in  equation  (20),  the  magnitude  (importance)  of  the  higher-order  terms, 
and  the  validity  of  the  assumption  that  p(k  and  gk]  are  independent  of  Ary  .  For 
sufficiently  small  Ary  ,  the  solutions  are  in  the  asymptotic  range  such  that  higher-order 
terms  are  negligible  and  the  assumption  that  p(k"  and  g[l>  are  independent  of  Ary.  is  valid. 

However,  achieving  the  asymptotic  range  for  practical  geometry  and  conditions  is  usually 
not  possible  and  m>3  is  undesirable  from  a  resources  point  of  view;  therefore,  methods  are 
needed  to  account  for  effects  of  higher-order  terms  for  practical  application  of  RE. 
Additionally,  methods  may  be  needed  to  account  for  possible  dependence  of  pyk  and  g[l> 

on  Ary  ,  although  not  addressed  herein.  Usually  5  *k  is  estimated  for  the  finest  value  of  the 
input  parameter,  i.e.,  8  *=8A*  corresponding  to  the  finest  solution  Sk  . 

For  777=3,  only  the  leading-order  term  can  be  evaluated.  Equations  are  obtained  for  SJ 
and  order-of-accuracy  pk 


'2k 


.Pk 


-1 


Pk 


ln(  8 


32u  lZ21t) 


ln(rk) 


(22) 

(23) 


Appendix  B  includes  results  for  777=5. 

Appendix  C  provides  verification  for  two  analytical  benchmarks  (one-dimensional 
wave  and  two-dimensional  Laplace  equations).  Multiple  solutions  were  used  to  evaluate 
the  RE  error  estimates,  including  the  effects  of  higher-order  terms.  Solving  for  the  first- 
order  term  is  relatively  easy  since  evaluation  of  equations  (22)  and  (23)  only  requires  that 
the  7?7=3  solutions  are  monotonically  convergent,  even  if  the  solutions  are  far  from  the 
asymptotic  range  and  equations  (22)  and  (23)  are  inaccurate.  Solving  for  the  higher-order 
terms  (i.e.,  second-order  term)  is  more  difficult  since  evaluation  of  the  777=5  solutions  for 
Sc,  Pk=h2) ,  and  gk=l2)  additionally  requires  that  the  solutions  are  relatively  close  to  the 
asymptotic  range,  i.e.,  within  about  6%  of  the  theoretical  order  of  accuracy  based  on  the 
modified  equation  pk  and  q k  . 
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The  solutions  show  that  equation  (22)  has  the  correct  form,  but  the  order  of 
accuracy  is  poorly  estimated  by  equation  (23)  except  in  the  asymptotic  range.  Therefore, 
one  approach  is  to  correct  equation  (22)  by  a  multiplication  correction  factor  to  account 
for  the  effects  of  higher-order  terms.  Two  correction  factors  were  investigated 


C 


k 


-1 

-1 


(24a) 


C 


k 


(e23t  /e 12,  -rqkk‘"){rkPk  -1)  |  (e23t  /e12t  -rkPt~)(rkp*  -1) 

|  j .  Pkest  ^  k est  ^  kk eSf  ^  Pkest  ^ .  ^kest  )^  ^kest  |  j 


(24b) 


pk  and  qk  are  estimates  for  the  1st  and  2nd  term  order  of  accuracy  pk '  and  /?{2) .  The 
estimated  values  can  be  based  either  on  pk^  and  q k  or  solutions  for  simplified  geometry 

and  conditions.  In  either  case,  preferably  including  the  effects  of  grid  stretching. 
Equation  (24a)  roughly  accounts  for  the  effects  of  higher-order  terms  by  replacing  pk  with 
pk  thereby  providing  an  improved  single-term  estimate.  Equation  (24b)  more 
rigorously  accounts  for  higher-order  terms  since  it  is  derived  from  the  two-term  estimate 
with  1st  and  2nd  term  order  of  accuracy  pk]  and  p[2)  replaced  by  pk  and  qk  .  Equation 
(24b)  simplifies  to  equation  (24a)  in  the  limit  of  the  asymptotic  range.  Both  correction 
factors  only  require  solutions  for  three  parameter  values.  Ck  <1  or  Ck  >  /  indicates  that 
the  leading-order  term  over  predicts  (higher-order  terms  net  negative)  or  under  predicts 
(higher-order  terms  net  positive)  the  error,  respectively.  Ck  given  by  equation  (24)  is 

fairly  universal  in  that  it  only  implicitly  depends  on  geometry  and  conditions.  However,  Ck 
is  based  on  results  from  only  two  linear  analytical  benchmarks  and  additional  benchmarks 
(especially  non-linear)  are  needed  to  confirm  the  universality  of  equation  (24)  or  to 
provide  alternative  forms. 


Combining  equation  (22)  and  (24)  provides  an  estimate  for  8  k  accounting  for  the 
effects  of  higher-order  terms 


Ct5’ 


C, 


'2K 


. Pk 


(25) 


The  estimate  includes  both  sign  and  magnitude.  Equation  (25)  is  used  to  estimate  U k  or 

8  *  and  Ukc  depending  on  how  close  the  solutions  are  to  the  asymptotic  range  (i.e.,  how 

close  Ck  is  to  1)  and  one’s  confidence  in  equation  (25).  There  are  many  reasons  for  lack 

of  confidence,  especially  for  complex  three-dimensional  flows.  Point  variables  invariably 
are  not  uniformly  convergent,  which  is  particularly  evident  near  inflection  points  and  zero 
crossings. 

Equations  (24)  and  (25)  need  further  testing  both  for  additional  analytical  benchmarks 
(as  already  mentioned)  and  practical  applications.  Also  alternative  strategies  for  including 
effects  of  higher-order  terms  may  be  just  as  viable.  Note  that  equation  (25)  differs 
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significantly  from  the  GCI  proposed  by  Roache  (1998).  Herein 
Ck  =  Ck  (e  ,rk,pk,pk  est ,  qk  est ) ,  whereas  in  the  GCI,  Ck  is  a  constant  referred  to  as  a 

factor  of  safety  Fs  which  equals  1.25  for  careful  grid  studies  and  3  for  cases  for  which  only 
two  grids  are  used. 

For  Ck  sufficiently  less  than  or  greater  than  1  and  lacking  confidence,  Uk  is 

estimated,  but  not  8  *k  and  Uk  .  Based  on  the  analytical  benchmark  studies  (Appendix  C), 

it  appears  that  equation  (25)  can  be  used  to  estimate  the  uncertainty  by  bounding  the  error 
by  the  sum  of  the  absolute  value  of  the  corrected  estimate  from  RE  and  the  absolute  value 
of  the  amount  of  the  correction 


Uk  = 


ck  5* 


+ 


(i-c*)s: 


(26) 


For  Ck  sufficiently  close  to  1  and  having  confidence,  8  *k  and  U k  are  estimated. 
Equation  (25)  is  used  to  estimate  the  error  8 *k ,  which  can  then  also  be  used  in  the 
calculation  of  Sc  [in  equation  (10)].  The  uncertainty  in  the  error  estimate  is  based  on  the 
amount  of  the  correction 


U,  = 


(1  -c,)S 


k  RE. 


(27) 


Note  that  in  the  limit  of  the  asymptotic  range,  Ck  =1,  8*  =  8*  =  ^*RE  ,  and  Uk  =0. 


3.2.4  Oscillatory  Convergence 

For  oscillatory  convergence,  i.e.,  condition  (ii)  in  equation  (16),  uncertainties  can 
be  estimated,  but  not  the  signs  and  magnitudes  of  the  errors.  Uncertainties  are  estimated 
based  on  determination  of  the  upper  ( SL/ )  and  lower  ( SL )  bounds  of  solution  oscillation, 

which  requires  more  than  m=3  solutions.  The  estimate  of  uncertainty  is  based  on  half  the 
solution  range 

(28) 


3.2.5  Divergence 

For  divergence,  i.e.,  condition  (iii)  in  equation  (16),  errors  and  or  uncertainties  can 
not  be  estimated.  The  preparation  and  verification  steps  must  be  reconsidered. 
Improvements  in  iterative  convergence,  parameter  specification  (e.g.,  grid  quality),  and/or 
CFD  code  may  be  required  to  achieve  converging  or  oscillatory  conditions. 


3.3  Validation 

Validation  is  defined  as  a  process  for  assessing  modeling  uncertainty  U SM  by  using 
benchmark  experimental  data  and,  when  conditions  permit,  estimating  the  sign  and 
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magnitude  of  the  modeling  error  SSM  itself.  Thus,  the  errors  and  uncertainties  in  the 

experimental  data  must  be  considered  in  addition  to  the  numerical  errors  and  uncertainties 
discussed  in  Section  3.2.  Approaches  to  estimating  experimental  uncertainties  are 
presented  and  discussed  by  Coleman  and  Steele  (1999). 

The  validation  methodology  of  Coleman  and  Stem  (1997)  which  properly  takes  into 
account  the  uncertainties  in  both  the  simulation  and  the  experimental  data  is  described  in 
this  section.  The  methodology  is  also  demonstrated  using  an  estimated  numerical  error 
and  corrected  simulation  and  validation  uncertainty  values. 


3.3.1  Methodology 

The  validation  comparison  for  a  simulated  and  measured  result  r  that  is  a  function 
of  the  variable  X  is  shown  in  figure  1.  The  experimentally  determined  r-value  of  the 
(A',.,?;.)  data  point  is  D  and,  as  before,  the  simulated  r-value  is  S.  Recall  from  equation  (1) 

that  the  simulation  error  8S  is  the  difference  between  S  and  the  truth  T.  Similarly,  the 
error  8D  in  the  data  is  the  difference  between  D  and  the  truth  T.  so  setting  the  simulation 
and  experimental  truths  equal  results  in 

D-Sd=S-8s  (29) 

The  comparison  error  E  is  defined  as  the  difference  of  D  and  S 

E  =  D- S  =d  D  -8S  =  8  D  -  (8 SMA  +Ss„  +S„)  (30) 

with  8sv/  decomposed  into  the  sum  of  ?>spd,  error  from  the  use  of  previous  data  such  as 
fluid  properties,  and  $sma,  error  from  modeling  assumptions.  Thus  E  is  the  resultant  of  all 
the  errors  associated  both  with  the  experimental  data  and  with  the  simulation.  For  the 
approach  in  which  no  estimate  8*SN  of  the  sign  and  magnitude  of  8  SN  is  made,  all  of  these 
errors  are  estimated  with  uncertainties.  (As  will  be  shown,  during  the  validation  process  an 
estimate  of  the  sign  and  magnitude  of  8  SMA  can  be  made  under  certain  conditions.) 


If  Xi,ri,  and  S  share  no  common  error  sources,  then  the  uncertainty  UE  in  the 
comparison  error  can  be  expressed  as 


Ul 


'})E 

2 

r dE ' 

u2d  + 

w 

u2=u2D+u2s 


(31) 


or 


u2  =u2 +u2  +u2  +u2 

u  E  u  D  ^  u  SMA  ~  u  SPD  '  u  SN 


(32) 


where  subscripts  are  used  in  the  same  manner  as  for  the  S's  . 

Ideally,  we  would  like  to  postulate  that  if  the  absolute  value  of  E  is  less  than  its 
uncertainty  UE ,  then  validation  is  achieved  (i.e.,  E  is  “zero”  considering  the  resolution 
imposed  by  the  “noise  level”  U E  ).  In  reality,  the  authors  know  of  no  approach  that  gives 
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an  estimate  of  U SMA ,  so  U E  cannot  be  estimated.  That  leaves  a  more  stringent  validation 
test  as  the  practical  alternative.  If  the  validation  uncertainty  Uv  is  defined  as  the 
combination  of  all  uncertainties  that  we  know  how  to  estimate  (i.e.,  all  but  USMA),  then 

U2  =U2 -U2  =u2 +u2  +u 2 

If  LEI  is  less  than  the  validation  uncertainty  Uv ,  the  combination  of  all  the  errors  in  D 
and  S  is  smaller  than  the  estimated  validation  uncertainty  and  validation  has  been  achieved 
at  the  Uv  level.  Uv  is  the  key  metric  in  the  validation  process.  Uv  is  the  validation 

“noise  level”  imposed  by  the  uncertainties  inherent  in  the  data,  the  numerical  solution,  and 
the  previous  experimental  data  used  in  the  simulation  model.  It  can  be  argued  that  one 
cannot  discriminate  once  LEI  is  less  than  this;  that  is,  as  long  as  LEI  is  less  than  this,  one 
cannot  evaluate  the  effectiveness  of  proposed  model  “improvements.” 

If  the  corrected  approach  of  equations  (3)-  (6)  is  used,  then  the  equations  equivalent 
to  equations  (30)  and  (33)  are 

Ec  —  D  —  Sc  —  8  D  — (8  SMA  +  8  SPD  +  £  SN )  (34) 


for  the  corrected  comparison  error  and 


-U 


SMA 


,TTZ  +TJ 1  +  TTZ 

D  ~  SPD  ~  SrN 


(35) 


for  the  corrected  validation  uncertainty.  Note  that  Sc  and  Ec  can  be  either  larger  or 
smaller  than  their  counterparts  S  and  E,  but  UE  and  Uv  should  be  smaller  than  U E  and 

Uv ,  respectively,  since  Us  N  should  be  smaller  than  USN  . 

For  the  data  point  (Xi ,ri),  UD  should  include  both  the  experimental  uncertainty  in  r 
and  the  additional  uncertainties  in  r  arising  from  experimental  uncertainties  in  the 
measurements  of  the  n  independent  variables  [x j )  in  X, .  The  expression  for  UD  that 
should  be  used  in  the  Uv  ( Uv  )  calculation  is  then 


Ul  =  U; 


+E 

j=i 


dr 

dx~. 

1  M 


(' ’J x, ): 


(36) 


In  some  cases,  the  terms  in  the  summation  in  equation  (36)  may  be  shown  to  be  very 
small,  using  an  order-of-magnitude  analysis,  and  then  neglected.  This  would  occur  in 
situations  in  which  the  Ux  values  are  of  "reasonable"  magnitude  and  gradients  in  r  are 

small.  In  regions  with  high  gradients  (e.g.,  near  a  surface  in  a  turbulent  flow),  these  terms 
may  be  very  significant  and  the  partial  derivatives  would  be  estimated  using  whatever 
{X.,r. )  data  is  available. 


There  is  also  a  very  real  possibility  that  measurements  of  different  variables  might 
share  identical  bias  errors.  This  is  easy  to  imagine  for  measurements  of  x,  y,  and  z. 
Another  possibility  is  D  and  S  sharing  an  identical  error  source,  for  example  if  the  same 
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density  table  (curve  fit)  is  used  both  in  data  reduction  in  the  experiment  and  in  the 
simulation.  In  such  cases,  additional  correlated  bias  terms  must  be  included  in  equation 
(31),  (32),  (33),  and  (35). 


To  estimate  U SPD  for  a  case  in  which  the  simulation  uses  previous  data  D  in  m 
instances,  one  would  need  to  evaluate 


U 


2 

SPD 


E 


'as  v 


dD, 


(37) 


where  the  UD  are  the  uncertainties  associated  with  the  data. 


3.3.2  Single  CFD  Code 

Consideration  of  equation  (32)  shows  that  (1)  the  more  uncertain  the  data,  and/or  (2) 
the  more  inaccurate  the  code  (greater  U SN  and  USPD),  the  easier  it  is  to  validate  a  code, 
since  the  greater  the  uncertainties  in  the  data  and  the  code  predictions,  the  greater  the 
noise  level  Uv .  However,  if  the  value  of  Uv  is  greater  than  that  designated  as  necessary 
in  a  research/design/development  program,  the  required  level  of  validation  could  not  be 
achieved  without  improvement  in  the  quality  of  the  data,  the  code,  or  both.  Also,  if  U SN 

and  U SPD  are  not  estimated,  but  I E\  is  less  than  U D ,  then  a  type  of  validation  can  be 

argued  to  have  been  achieved,  but  clearly  as  shown  by  the  present  methodology,  at  an 
unknown  level. 

If  there  is  a  programmatic  validation  requirement,  denote  it  as  U  d  since  validation  is 
required  at  that  uncertainty  level  or  below.  From  a  general  perspective,  if  we  consider  the 
three  variables  Uv ,  \E\ ,  and  U,  there  are  six  combinations  (assuming  none  of  the  three 
variables  are  equal): 

1.  \E\<  Uv  <  Ureqd 

2.  \E\  <  Ureqd  <  Uv 

3.  Ureqd  <  \E\  <  Uv 

4.  Uy  <  \E\  <  Ureqd  (38) 

5-  Uv<  Ureqd  <  \E\ 

6-  Ureqd  <UV<  \E\ 

In  cases  1,  2  and  3,  \E\<UV;  validation  is  achieved  at  the  Uv  level;  and  the 
comparison  error  is  below  the  noise  level,  so  attempting  to  decrease  the  error  6  SMA  due  to 
the  modeling  assumptions  in  the  simulation  is  not  feasible  from  an  uncertainty  standpoint. 
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In  case  1,  validation  has  been  achieved  at  a  level  below  U  d,  so  validation  is  successful 
from  a  programmatic  standpoint. 

In  cases  4,  5  and  6,  t/K<|is| ,  so  the  comparison  error  is  above  the  noise  level  and 
using  the  sign  and  magnitude  of  E  to  estimate  8  SMA  is  feasible  from  an  uncertainty 
standpoint.  If  UV«\E\,  then  E  corresponds  to  8 SMA  and  the  error  from  the  modeling 
assumptions  can  be  determined  unambiguously.  In  case  4,  validation  is  successful  at  the 
\E\  level  from  a  programmatic  standpoint. 

A  similar  comparison  table  can  be  constructed  using  \EC  I,  Uv  ,  and  Ureqd.  Since  Ec 
can  be  larger  or  smaller  than  E,  but  Uv  should  always  be  less  than  Uv ,  the  results  for  a 

given  corrected  case  are  not  necessarily  analogous  to  those  for  the  corresponding 
uncorrected  case.  That  is,  a  variable  can  be  validated  in  the  corrected  but  not  in  the 
uncorrected  case,  or  vice  versa.  However,  the  band  Ec  ±  UE  should  always  give  a 

smaller  (therefore  better)  range  within  which  the  true  value  of  E  lies  than  the  band  E  ±  Ue, 
assuming  that  one’s  confidence  in  using  the  estimate  d*SN  is  not  misplaced.  Furthermore, 
for  cases  4,  5,  and  6,  one  can  argue  that  Ec  more  likely  corresponds  to  8  SMA . 

In  general,  validation  of  a  code's  predictions  of  a  number  (N)  of  different  variables  is 
desired,  and  this  means  that  in  a  particular  validation  effort  there  could  be  N  different  E, 
Ec,  Uv ,  Uv  ,  and  Ureqd  values  and  (perhaps)  some  successful  and  some  unsuccessful 
validations.  For  each  variable,  a  plot  of  the  simulation  prediction  versus  X  compared  with 
the  {Xt ,  r )  data  points  gives  a  traditional  overview  of  the  validation  status,  but  the 
interpretation  of  the  comparison  is  greatly  affected  by  choice  of  the  scale  and  the  size  of 
the  symbols.  A  plot  of  ±UV(±  U  Vr  )  and  E  (Ec),  and  Ureqd  (if  known)  versus  X  for  each 

variable  is  particularly  useful  in  drawing  conclusions,  and  the  interpretation  of  the 
comparison  is  more  insensitive  to  scale  and  symbol  size  choices. 


3.3.3  Comparison  of  Multiple  Codes  and/or  Models 

When  a  validation  effort  involves  multiple  codes  and/or  models,  the  procedure 
discussed  above  —  comparison  of  values  of  E  and  Uv  (and  U  d  if  known)  for  the  N 
variables  —  should  be  performed  for  each  code/model. 

Since  each  code/model  may  have  a  different  Uv ,  some  method  to  compare  the 
different  codes’/models’  performance  for  each  variable  in  the  validation  is  useful.  The 
range  within  which  (95  times  out  of  100)  the  true  value  of  E  lies  is  E±Ue.  From 
equation  (32),  when  USMA  is  zero  then  Uv  =  UE  ,  so  for  that  ideal  condition  the  maximum 
absolute  magnitude  of  the  95%  confidence  interval  is  given  by  \E\  +  Ur.  Comparison  of 
the  (\e\  +  Uv  Y s  for  the  different  codes/models  then  shows  which  has  the  smallest  range  of 
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likely  error  assuming  all  USMA’ s  are  zero.  This  allows  appropriate  comparisons  of  (low 
is)/(high  UY )  with  (high  E/low  Uv )  codes/modcls. 

A  similar  discussion  holds  if  the  corrected  values  are  used. 


3.3.4  Predictions  of  Trends 

In  some  instances,  the  ability  of  a  code  or  model  to  predict  the  trend  of  a  variable  may 
be  the  subject  of  a  validation  effort.  An  example  would  be  the  difference  in  drag  for  two 
ship  configurations  tested  at  the  same  Froude  number.  The  procedure  discussed  above  — 
comparison  of  LEI  and  Uv  for  the  drag  —  should  be  performed  for  each  configuration.  The 

difference  A  in  drag  for  the  two  configurations  should  then  be  considered  as  the  variable 
that  is  the  subject  of  the  validation.  As  discussed  in  Coleman  and  Steele  (1999),  because 
of  correlated  bias  uncertainty  effects  in  the  experimental  data  the  magnitude  of  the 
uncertainty  in  A  may  be  significantly  less  than  the  uncertainty  in  either  of  the  two 
experimentally  determined  drag  values.  This  means  that  the  value  of  Uv  for  A  may  be 

significantly  less  than  the  Uv 's  for  the  drag  values,  allowing  for  a  more  stringent  validation 

criterion  for  the  difference  than  for  the  absolute  magnitudes  of  the  variables.  Choice  of  the 
corrected  or  uncorrected  approach  should  be  made  on  a  specific  case-by-case  basis. 


3.3.5  Corrected  vs.  Uncorrected  Simulation  Results 

If  a  validation  using  the  corrected  approach  is  successful  at  a  set  condition,  then  if  one 
chooses  to  associate  that  validation  uncertainty  level  with  the  simulation's  prediction  at  a 
neighboring  condition  that  prediction  must  also  be  corrected.  That  means  enough  runs  are 
required  at  the  new  condition  to  allow  estimation  of  the  numerical  errors  and 
uncertainties.  If  this  is  not  done,  then  the  comparison  error  E  and  validation  uncertainty 
Ur  corresponding  to  the  use  of  the  uncorrected  S  and  its  associated  (larger)  USN  should 

be  the  ones  considered  in  the  validation  with  which  one  wants  to  associate  the  prediction 
at  a  new  condition.  (Whether  to  and  how  to  associate  an  uncertainty  level  at  a  validated 
condition  with  a  prediction  at  a  neighboring  condition  is  very  much  unresolved  and  is 
justifiably  the  subject  of  much  debate  at  this  time.) 

As  discussed  in  Section  3.3.2,  however,  the  band  Ec  ±  UE  should  always  give  a 

smaller  (therefore  better)  range  within  which  the  true  value  of  E  lies  than  the  band  E  ±  Ue, 
assuming  that  one’s  confidence  in  using  the  estimate  6  *^  is  not  misplaced. 


17 


4.0  Example  for  RANS  CFD  Code 


Example  results  of  verification  and  validation  are  presented  for  a  single  CFD  code  and 
for  specified  objectives,  geometry,  conditions,  and  available  benchmark  information.  The 
CFD  code  is  CFDSHIP-IOWA,  which  is  a  general-purpose,  multi-block,  high  performance 
computing  (parallel),  unsteady  RANS  code  (Paterson  et  al,  1998;  Wilson  et  al.,  1998) 
developed  for  computational  ship  hydrodynamics.  The  RANS  equations  are  solved  using 
higher-order  upwind  finite  differences,  PISO,  Baldwin-Fomax  turbulence  model,  and 
exact  and  approximate  treatments,  respectively,  of  the  kinematic  and  dynamic  free-surface 
boundary  conditions.  The  objectives  are  to  demonstrate  the  usefulness  of  the  proposed 
verification  and  validation  procedures  and  methodology  and  establish  the  levels  of 
verification  and  validation  of  the  simulation  results  for  an  established  benchmark  for  ship 
hydrodynamics  CFD  validation. 


4.1  Geometry,  Conditions,  and  Benchmark  Data 

The  geometry  is  the  Series  60  cargo/container  ship.  The  Series  60  was  used  for  two 
of  the  three  test  cases  at  the  last  international  workshop  on  validation  of  ship 
hydrodynamics  CFD  codes  (CFD  Workshop  Tokyo,  1994).  The  conditions  for  the 
calculations  are  Froude  number  Fr  =  0.316,  Reynolds  number  Re  =  4.3xl06,  and  zero 
sinkage  and  trim.  These  are  the  same  conditions  as  the  experiments,  except  the  resistance 
and  sinkage  and  trim  tests,  as  explained  next.  The  variables  selected  for  verification  and 
validation  are  resistance  CT  (integral  variable)  and  wave  profile  'Q  (point  variable). 

The  benchmark  data  is  provided  by  Toda  et  al.  (1992),  which  was  also  the  data  used 
for  the  Series  60  test  cases  at  the  CFD  Workshop  Tokyo  (1994).  The  data  includes 
resistance  and  sinkage  and  trim  for  a  range  of  Fr  for  the  model  free  condition  (i.e.,  free  to 
sink  and  trim);  and  wave  profiles,  near-field  wave  pattern,  and  mean  velocities  and 
pressures  at  numerous  stations  from  the  bow  to  the  stem  and  near  wake,  all  for  Fr  = 
(0.16,  0.316)  and  the  zero  sinkage  and  trim  model  fixed  condition.  The  data  also  includes 
uncertainty  estimates,  which  were  recently  confirmed/updated  by  Fongo  and  Stern  (1999) 
closely  following  standard  procedures  (Coleman  and  Steele,  1999). 

The  resistance  is  known  to  be  larger  for  free  vs.  fixed  models.  Data  for  the  Series 
60  indicates  about  an  8%  increase  in  Ct  for  the  free  vs.  fixed  condition  over  a  range  of  Fr 
including  F/ -0.3 1 6  (Ogiwara  and  Kajatani,  1994).  The  Toda  et  al.  (1992)  resistance 
values  were  calibrated  (i.e.,  reduced  by  8%)  for  effects  of  sinkage  and  trim  for  the  present 
comparisons. 

4.2  Computational  Grids 

Grid  studies  were  conducted  using  four  grids  (m= 4),  which  enables  two  separate 
grid  studies  to  be  performed  and  compared.  Grid  study  1  gives  estimates  for  grid  errors 
and  uncertainties  on  grid  1  using  the  three  finest  grids  1-3  while  grid  study  2  gives 
estimates  for  grid  errors  and  uncertainties  on  grid  2  using  the  three  coarsest  grids  2-4.  The 
results  for  grid  study  1  are  given  in  detail  and  the  differences  for  grid  study  2  are  also 
mentioned.  The  grids  were  generated  using  the  commercial  code  GRIDGEN  (Pointwise, 
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Inc.)  with  consideration  to  topology;  number  of  points  and  grid  refinement  ratio  rG;  near¬ 
wall  spacing  and  turbulence  model  requirement  that  first  point  should  be  at  y+<  1 ;  bow  and 
stem  spacing;  and  free-surface  spacing. 

The  topology  is  body-fitted,  H- type,  and  single  block..  The  sizes  of  grids  1  (finest) 
through  4  (coarsest)  are  287x78x43  =  876,211,  201x51x31  =  317,781,  144x36x22  = 

114,048,  and  101x26x16  =  42,016,  and  the  grid  refinement  ratio  rG  =  42 .  Clustering  was 

used  near  the  bow  and  stern  in  the  ^-direction,  at  the  hull  in  the  T| -direction,  and  near  the 
free  surface  in  the  ^-direction.  The  v+  values  for  grids  1-4  were  about  0.7,  1,  1.4,  and  2, 
respectively.  About  twice  the  number  of  grid  points  in  the  r| -direction  would  be  required 
to  achieve  v+  <1.0  for  all  four  grids  1-4  (i.e.,  roughly  1,800,000  points  on  the  finest  grid). 

With  grid  refinement  ratio  rG  =  42 ,  only  grids  1  and  2  were  generated.  Grids  3  and  4 

were  obtained  by  removing  every  other  point  from  grids  1  and  2,  respectively  (i.e.,  the 
grid  spacing  of  grids  3  and  4  is  twice  that  of  grids  1  and  2,  respectively).  Grids  1  and  2 
were  generated  by  specifying  the  grid  spacing  at  the  corners  and  number  of  points  along 
the  edges  of  the  computational  blocks.  The  faces  of  the  computational  blocks  were 
smoothed  using  an  elliptic  solver  after  which  the  coordinates  in  the  interior  were  obtained 
using  transfinite  interpolation  from  the  block  faces.  Grid  2  was  generated  from  grid  1  by 
increasing  the  grid  spacing  and  decreasing  the  number  of  computational  cells  in  each 
coordinate  direction  at  the  corners  of  the  blocks  by  a  factor  rG.  A  comparison  of  the  four 
grids  at  the  free  surface  plane  is  shown  in  figure  2  along  with  computed  wave  elevation 
contours 


4.3  Verification  and  Validation  of  Integral  Variable:  Resistance 

Verification.  Verification  was  performed  with  consideration  to  iterative  and  grid 
convergence  studies,  i.e.,  8 SN  =  S;  +SG  and  U2SN  =  U2r  +UG. 

Iterative  convergence  was  assessed  by  examining  iterative  history  of  ship  forces  and 
L2  norm  of  solution  changes  summed  over  all  grid  points.  Figure  3  shows  a  portion  of  the 
iterative  history  on  grid  1.  The  portion  shown  represents  a  computation  started  from  a 
previous  solution  and  does  not  reflect  the  total  iterative  history.  Solution  change  drops 
four  orders  of  magnitude  from  an  initial  value  of  about  10 2  (not  shown)  to  a  final  value  of 
10"6.  The  variation  in  CY  is  about  0.07 %SG  over  the  last  period  of  oscillation  (i.e.,  Ui  = 
0. 07 %SG).  Iterative  uncertainty  is  estimated  as  half  the  range  of  the  maximum  and 
minimum  values  over  the  last  two  periods  of  oscillation  (see  figure  3c).  Iterative  histories 
for  grids  2-4  show  iterative  uncertainties  of  about  0.02,  0.03,  and  0.01  %SG,  respectively. 
The  level  of  iterative  uncertainties  Ui  for  grids  2-4  are  at  least  two  orders  of  magnitude 
less  than  the  corresponding  grid  uncertainties  UG,  whereas  the  iterative  uncertainty  for 
grid  1  is  only  one  order  of  magnitude  smaller  than  the  grid  error.  For  all  four  grids  the 
iteration  errors  and  uncertainties  are  assumed  to  be  negligible  in  comparison  to  the  grid 
errors  and  uncertainties  for  all  four  solutions  (i.e.,  8/  «  SG  and  Ui«  UG  such  that  8s, v  = 
SG  and  Usn=Ug ). 

The  results  from  the  grid  convergence  study  for  CY  are  summarized  in  tables  2  and  3. 
The  solutions  for  CY  indicate  the  converging  condition  (i)  of  equation  (16)  with 
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Rg  =e2l  /e,,  =0.21 .  The  first-order  RE  estimate  8 REc  [in  equation  (22)],  order  of  accuracy 
pG  [in  equation  (23)],  and  correction  factor  CG  [in  equation  (24a)]  are 


RErn 


(  g  3 

fc21G 

'  0.06x1 0‘3  "j 

rPa  -1 

VG  J 

[(72)“ -1 J 

0.02x10 


-3 


Pg 


ln(£32G/82lG)_ln(0.28/0.06) 


ln(rG) 


ln(>/2) 


:  4.4 


(39) 

(40) 


,  _  C  - 1  _(V2)44-1_ 

G  r^'-l  (a/2)2-1 


(41) 


where  pest=pth= 2  was  used  in  equation  (41).  Uncertainty  and  error  estimates  are  made 
next  both  considering  CG  as  sufficiently  less  than  or  greater  than  1  and  lacking  confidence 
and  CG  as  close  to  1  and  having  confidence,  as  discussed  in  Section  3.2.3. 

For  CG  =  3.7  considered  as  sufficiently  less  than  or  greater  than  1  and  lacking 
confidence,  UG  is  estimated  and  not  8C 

Ua  =  |  + 1(1  —  CG )5  ;Foi  |  =  0.06.y1  0'3  +  0.05 xl  0~3  =  0. 1  lxl  O"3  (42) 


UG  is  2. 1  %  SC| . 


For  CG  =  3.7  considered  close  to  1  and  having  confidence,  both  and  8*  and  UGc  are 
estimated 


S;  =Cg8^g  =  0.06x1 0“3 


(43) 


uo,  = 


(i-cG)s: 


=  0.05xl0~ 


(44) 


The  corrected  solution  SG  is  defined  with  S  =  SGi 

Sc  =SGi  -8Gi  =4.99x10-3  (45) 

8  G|  and  UGc  are  1.2%  and  1.0%  Sc,  respectively.  In  both  cases,  the  level  of  verification  is 
relatively  small  <2. 1  %  5*Gi  . 

Table  3  includes  results  for  grid  study  2,  which  are  similar  to  those  for  grid  study  1, 
but  the  values  are  larger  by  a  factor  of  about  3,  except  SG  which  differs  by  only  3%.  Also 
shown  in  table  2  are  the  pressure  Cp  and  frictional  CV  components  of  Ct ■  Cf  comprises 
about  70%  of  CT  and  also  displays  convergence;  however,  CP  is  convergent  for  the  second 
grid  study  and  neutrally  convergent  (RG= 0)  to  three  significant  figures  for  the  first  grid 
study  (i.e.,  Cp  is  grid  independent  on  the  finest  grid).  Solution  changes  between  grids  1 
and  2  for  Cp  are  at  or  below  the  level  of  iteration  uncertainty  (0.1%S'G),  so  that  further 

grid  refinement  is  unwarranted.  Apparently  for  this  geometry,  convergence  of  CV  with  grid 
refinement  is  slower  than  that  of  CP.  The  results  show  that  the  use  of  finer  grids  is 

problematic;  since,  the  next  largest  grid  with  rG  =  V2  would  have  2.4M  grid  points  and 
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iterative  errors  and  grid  errors  would  likely  be  of  similar  order  of  magnitude.  The  RG,  pG, 
and  CG  values  are  far  from  their  asymptotic  range  values  of  RG=0.5,  pG= 2,  and  CG~  1, 
respectively. 


Validation.  Validation  is  performed  using  both  the  simulation  prediction  S  and  the 
corrected  simulation  prediction  SG,  as  summarized  in  table  4.  First  using  S,  the  comparison 
error  is  calculated  from  equation  (30)  with  S  =  SG  as 

E  =  D-S  =  5 .42x1  (T3  - 5.05xl(T3  =  0.37x10  3  =  6.8 %D  (46) 

The  validation  uncertainty  is  calculated  from  equation  (33)  as 

Uv  =^U2SN+Ul  =  0.17x10“3  =3.1%£)  (47) 

where  Usn=Ug  =1.9 %D  and  Uo=2.5%D.  Comparison  error  \E\  >Uv  such  that  the 

simulation  results  are  not  validated.  Usn  and  Ud  are  of  similar  order  such  that  reduction  in 
Uv  would  require  reduction  of  UD  and  Usn-  Reduction  of  Usn  by  using  finer  grids  may  be 
possible;  however,  as  already  mentioned,  iterative  errors  will  likely  be  of  similar  order  of 
magnitude  and  will  also  need  to  be  accurately  estimated.  E  is  positive,  i.e.,  the  simulation 
under  predicts  the  data.  The  trends  shown  in  table  2  suggest  Cp  too  small.  Presumably 
modeling  errors  such  as  resolution  of  the  wave  field  and  inclusion  of  effects  of  sinkage  and 
trim  can  be  addressed  to  reduce  E  and  validate  CT  at  TV=3. 1  %D;  however,  the  case  for 
this  reasoning  is  stronger  when  considering  the  corrected  comparison  error,  as  discussed 
next. 

Second  using  Sc,  the  corrected  comparison  error  is  calculated  from  equation  (34)  as 

Ec  =  £>-Sc  =5.42x10  3  -4.99xl(T3  =  0.43xl0~3  =  7.9%£>  (48) 

The  validation  uncertainty  is  calculated  from  equation  (35)  as 

UVc  =  ^UlcN+Ul  =  0. 14x1 0'3  =  2.6  %D  (49) 

where  Us .  N  =  UG  =0.8% 74  Here  again,  | EC\>UV  such  that  the  simulation  results  are 
not  validated.  However,  validation  uncertainty  Uv  is  relatively  small  and  Us  N«UD 
more  strongly  suggests  than  was  the  case  for  E  that  Ec  is  mostly  due  to  modeling  errors. 
Therefore  modeling  issues  should/can  be  improved  to  reduce  Ec  and  validate  Ct  at  the 
reduced  level  Uv  =2.6 %D  in  comparison  to  equation  (47). 

The  results  from  grid  study  2  are  summarized  in  table  5.  Note  that  validation  of  the 
comparison  error  E  is  achieved  at  the  level  of  Uv=6.7%D  while  validation  of  the  corrected 
comparison  error  Ec  is  not. 


4.4  Verification  and  Validation  of  a  Point  Variable:  Wave  Profile 

Verification.  Verification  for  the  wave  profile  was  conducted  as  per  that 
described  for  the  resistance  in  Section  4.3  with  the  distinction  that  a  point  variable  is 
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defined  over  a  distribution  of  grid  points.  Interpolation  of  the  wave  profile  on  all  grids 
onto  a  common  distribution  is  required  to  compute  solution  changes.  Since  calculation  of 
the  comparison  error  E=D-S  is  required  for  validation,  wave  profiles  on  grids  1-4  are 
interpolated  onto  the  distribution  of  the  data.  The  same  four  grids  were  used  and,  here 
again  iteration  errors  and  uncertainties  were  negligible  in  comparison  to  the  grid  errors 
and  uncertainties  for  all  four  solutions,  i.e.,  8/  «  Sc  and  Uj«  UG  such  that  8SN  =  Sc  and 
Usn=Ug. 


Rg  at  local  maximums  and  minimums  (i.e.,  x/L  =  0.1,  0.4,  and  0.65  in  figure  4a) 
and  based  on  L2  norm  solution  changes  both  show  convergence.  The  spatial  order  of 
accuracy  for  the  wave  profile  was  computed  from  the  L2  norm  of  solution  changes 


In  (rG) 


=  1.3 


(50) 


where  <  >  is  used  to  denote  a  profile-averaged  value  and  ||e  ||0  denotes  the  L2  norm  of 
solution  change  over  the  N  points  in  the  region,  0  <  x/L  <  1 


(51) 


Correction  factor  is  computed  from  equation  (24a)  using  order  of  accuracy  pG  in  equation 
(50)  and  pc  =  2.0 

^est 


Jpg)  _i 

G _ / 

yPGest  _1 

'G  1 


(V2)u-1 

(V2)2-l 


0.56 


(52) 


The  estimates  for  order  of  accuracy  and  correction  factor  in  equations  (50)  and  (51)  were 
used  to  estimate  grid  error  and  uncertainty  for  the  wave  profile  at  each  grid  point. 

For  <CG>  =  0.56  considered  as  sufficiently  less  than  or  greater  than  1  and  lacking 
confidence,  pointwise  values  for  UG  are  estimated  and  not  bG.  Equation  (26)  is  used  to 
estimate  UG 


Ur 


(CG 


■'21,. 


Ipg) 


-1 
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«-{C0>) 


'21,- 


Ipg) 


-1 


(53) 


For  <Cg>=0.56  considered  close  to  1  and  having  confidence,  pointwise  values  for  both 
8  *  and  UG  are  estimated  using  equations  (25)  and  (27) 


K=(Cg 


'21, 


Apg) 


Ua  = 


a  -<Q>) 


'21,. 


■,(Pg) 


(54) 

(55) 


Equation  (10)  is  used  to  calculate  Sc  at  each  grid  point 
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(56) 


5 


c 


The  results  are  summarized  in  table  6.  The  level  of  verification  is  similar  to  that 
for  Ct  with  slightly  higher  values.  Table  6  includes  results  for  grid  study  2,  which  are 
closer  to  those  for  grid  study  1  than  was  the  case  for  Ct,  i.e.,  are  only  larger  by  a  factor  of 
2  vs.  3  for  CT.  The  RG,  pG,  and  CG  values  are  closer  to  and  seem  to  be  approaching  the 
asymptotic  range. 


Validation.  Validation  of  the  wave  profile  is  performed  using  both  the  simulation 
prediction  S  and  the  corrected  simulation  prediction  S(  .  Profile- averaged  values  for  both 
definitions  of  the  comparison  error,  validation  uncertainty,  and  simulation  uncertainty  are 
given  in  table  7.  Values  are  normalized  with  the  maximum  value  for  the  wave  profile 
Cmax= 0.014  and  the  uncertainty  in  the  data  was  reported  to  be  3.7 %C,max.  For  grid  study  1, 
E  is  nearly  validated  at  about  5%.  The  trends  are  similar  to  those  for  Ct,  except  there  are 
smaller  differences  between  the  use  of  E  and  Ec. 

The  point  comparison  error  E=D-S  is  compared  to  validation  uncertainty  Uy  in 
figure  4b,  while  error  EC=D-SG  is  compared  to  validation  uncertainty  Uy  in  figure  4d.  In 
the  latter  case,  the  validation  uncertainty  Uy  in  figure  4d  is  mostly  due  to  U[>.  Much  of  the 
profile  is  validated.  The  largest  errors  are  at  the  crests  and  trough  regions,  i.e.,  bow, 
shoulder,  and  stem  waves. 

The  results  from  grid  study  2  are  summarized  in  table  8  and  included  in  Figure  4. 
The  results  are  similar  to  those  for  grid  study  1,  but  both  E  and  Ec  and  Uy  and  Uv<  are 
larger. 


5.  Conclusions  and  Recommendations 

The  verification  and  validation  procedures  and  methodology  presented  should  have 
applicability  to  a  fairly  broad  range  of  CFD  codes,  including  RANS,  Navier- Stokes,  Euler, 
boundary-element  methods,  and  others.  The  concepts  and  definitions  and  associated 
mathematical  framework  are  well  founded.  However,  clearly  much  more  work  is  needed 
for  other  CFD  codes  (such  as  large-eddy  simulations),  additional  error  sources,  and 
alternative  error  and  uncertainty  estimation  methods,  e.g.,  single-grid  methods  and 
alternative  strategies  to  account  for  the  effects  of  higher-order  terms  in  RE.  Furthermore, 
more  experience  is  needed  through  application  for  different  codes  and  geometry  and 
conditions. 

Nonetheless,  the  verification  and  validation  procedures  and  methodology  are 
recommended  for  use.  Use  of  such  procedures  and  methodology  should  be  helpful  in 
guiding  future  developments  in  CFD  through  documentation,  verification,  and  validation 
studies  and  in  transition  of  CFD  codes  to  design  through  establishment  of  credibility. 
Presumably,  with  a  sufficient  number  of  documented,  verified,  and  validated  solutions 
along  with  selected  verification  studies  a  CFD  code  can  be  accredited  for  a  certain  range 
of  applications.  The  contribution  of  the  present  work  is  in  providing  procedures  and 
methodology  for  the  former,  which  hopefully  will  help  lead  to  the  latter. 
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Table  1.  Definitions  of  errors  and  uncertainties  and  verification  and  validation 


Errors 

Present  and  Roache 
(1998) 

Error  8  is  the  difference  between  a  simulation  value  or  an 
experimental  value  and  the  truth 

AIAA  (1998) 

A  recognizable  deficiency  in  any  phase  or  activity  of  modeling  and 
simulation  that  is  not  due  to  lack  of  knowledge 

Uncertainties 

Present  and  Roache 
(1998) 

An  uncertainty  U  is  an  estimate  of  an  error  such  that  the  interval 
±U  contains  the  true  value  of  8  95  times  out  of  100 

AIAA  (1998) 

A  potential  deficiency  in  any  phase  or  activity  of  the  modeling 
process  that  is  due  to  lack  of  knowledge 

Verification 

Present 

Verification  is  defined  as  a  process  for  assessing  numerical 
uncertainty  USN  and,  when  conditions  permit,  estimating  the  sign 

and  magnitude  of  the  numerical  error  8  *SN  itself  and  the  uncertainty 
Us  N  in  that  error  estimate. 

Roache (1998) 

Solving  the  equations  right/mathematics 

AIAA  (1998) 

The  process  of  determining  that  a  model  implementation  accurately 
represents  the  developer’s  conceptual  description  of  the  model  and 
the  solution  to  the  model 

Validation 

Present 

Validation  is  defined  as  a  process  for  assessing  modeling 
uncertainty  USM  by  using  benchmark  experimental  data  and,  when 
conditions  permit,  estimating  the  sign  and  magnitude  of  the 
modeling  error  8  SM  itself. 

Roache (1998) 

Solving  the  right  equations/science/engineering 

AIAA  (1998) 

The  process  of  determining  the  degree  to  which  a  model  is  an 
accurate  representation  of  the  real  world  from  the  perspective  of 
the  intended  uses  of  the  model 
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Table  2  Grid  convergence  study  for  total  CT,  pressure  CP,  and  frictional  CF 
resistance  (xlO3)  for  Series  60. 


Grid 

Grid  4 

Grid  3 

Grid  2 

Grid  1 

Data 

101x26x16 

144x36x22 

201x51x31 

287x71x43 

CT 

6.02 

5.39 

5.11 

5.05 

5.42 

£ 

-10% 

-5.2% 

-1.2% 

Cp 

1.88 

1.60 

1.60 

CR  =  2.00 

£ 

-0.6% 

0.0% 

cF 

4.14 

3.69 

3.51 

3.45 

3.42 

£ 

-11% 

-4.9% 

-1.7% 

ITTC 

%SG. 


Table  3.  Verification  of  total  resistance  CF  (xlO  3)  for  Series  60. 


Study 

Rg 

Pg 

cG 

uG 

K 

Ur 

GC 

Sc 

1 

(grids  1-3) 

0.21 

4.4 

3.7 

2.1% 

1.2% 

0.9% 

4.99 

2 

(grids  2-4) 

0.44 

2.3 

1.3 

6.7% 

5.5% 

1.1% 

4.83 

%SG. 
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Table  4.  Validation  of  total  resistance  for  Series  60  -  study  1  (grids  1-3). 


E% 

Uv% 

UD% 

Usn% 

E=D-S 

6.8 

3.1 

2.5 

1.9 

Ec=D-Sc 

7.9 

2.6 

2.5 

0.8 

%D. 


Table  5.  Validation  of  total  resistance  for  Series  60  -  study  2  (grids  2-4). 


E% 

Uv% 

UD% 

UsN% 

E=D-S 

5.7 

6.7 

2.5 

6.3 

Ec=D-Sc 

11 

2.7 

2.5 

1.0 

%D. 


Table  6  Profile-averaged  values  from  verification  of  wave  profile  for  Series  60. 


Study 

Rg 

Pg 

cG 

uG 

Uc 

GC 

1 

(grids  1-3) 

0.64 

1.3 

0.56 

2.0% 

0.9% 

2 

(grids  2-4) 

0.68 

1.1 

0.47 

4.1% 

2.2% 

>max  • 


Table  7.  Profile-averaged  values  from  validation 
of  wave  profile  for  Series  60  -  study  1  (grids  1-3). 


E% 

Uv% 

UD% 

UsN% 

E=D-S 

5.2 

4.2 

3.7 

2.0 

Ec=D-Sc 

5.6 

3.8 

3.7 

0.9 

o/r 

/  O^rnax  * 


Table  8.  Profile-averaged  values  from  validation 
of  wave  profile  for  Series  60  -  study  2  (grids  2-4). 


E% 

Uv% 

UD% 

Usn% 

E=D-S 

5.6 

5.5 

3.7 

4.1 

Ec=D-Sc 

6.6 

4.3 

3.7 

2.2 

28 
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Figure  1  Definition  of  comparison  error. 


29 


Figure  2.  Grids  and  wave  contours  from  verification  and  validation  studies  for  Series  60: 

(a)  and  (b)  coarsest  -  grid  4;  (c)  and  (d)  grid  3;  (e)  and  (f)  grid  2;  and  (g)  and  (h) 
finest  -  grid  1. 


30 


Iteration 


Figure  3.  Iteration  history  for  Series  60  on  grid  1:  (a)  solution  change,  (b)  ship  forces  -  CF, 
CP,  and  CT  and  (c)  magnified  view  of  total  resistance  CT  over  last  two  periods  of 
oscillation. 
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Appendix  A.  Derivation  of  Simulation  Error  Equation 

There  are  three  engineering  approaches  to  solving  fluid  mechanics  problems: 
analytical,  experimental,  and  simulation.  Fluid  mechanics  problems  are  governed  by  the 
laws  of  physics,  which  are  formulated  for  unsteady  flow  as  initial  and  boundary  value 
problems  (IB VP)  and  for  steady  flow  problems  as  boundary  value  problems  (BVP). 

The  IB  VP  is  defined  by  a  continuous  partial  differential  equation  (PDE)  operator 
Lt  with  specified  initial  (IC)  and  boundary  (BC)  conditions 

Lt(T)=  0 

IC  :  T(x,t  =  0)  =  Gt (x)  (A.l) 

BC  :  T (xB ,  t)  =  Ht  (t) 

x  is  the  spatial  coordinate(s)  and  may  be  a  vector,  the  functions  Gr  and  IIT  are  the  IC  (at 
t=0)  and  BC  (at  x  =  xB),  respectively,  t  is  time,  and  T  is  the  true  or  exact  solution.  By 
definition,  equation  (A.l)  contains  no  modeling  or  numerical  errors. 

The  experimental  approach  does  not  solve  equation  (A.l),  but  instead  uses 
experimental  measurement  systems  to  determine  T.  This  process  results  in  bias  and 
precision  errors  that  lead  to  an  uncertainty  Ud  in  the  experimental  measurement  D. 

Analytical  and  simulation  approaches  formulate  the  IB  VP  by  selection  of  the  PDE, 
IC,  and  BC  to  model  the  physical  phenomena 

Lm(M)  =  0 

IC  :  M(x,  t  =  0)  =  Gm  (jc)  (A.  2) 

BC :  M(xB,t)  =  HM  (t) 

with  similar  definitions  as  per  equation  (A.l);  however,  LM,  GM,  Hm,  and  xB  all  may 
contain  modeling  assumptions  such  that  M  . 

Assumptions  are  made  in  modeling  geometry,  turbulence,  non-Newtonian  fluids, 
combustion,  compressibility,  two-fluid  and  rarified  gas  flows,  etc.  An  IBVP  for  the 
modeling  error  8  SM  =M  -T  (i.e.,  the  difference  between  the  model  and  true  values)  can 

be  obtained  by  subtracting  equation  (A.l)  and  (A. 2),  then  subtracting  LM  (T)  from  both 
sides  of  that  result,  and  lastly  assuming  that  the  operator  LM  is  linear 

lm  C M-T)  =  lm  (8  sm  )  =  rM=  -lm  C t ) 
ic :  8  SM  (x,  t  =0)  =  Gm  (x)  -  Gt  (x)  (A.3) 

BC:8SM(xB,t)  =  HM(t)-HT(t) 

The  assumption  that  LM  is  linear  is  a  major  limitation  since  most  fluid  mechanics  problems 
of  interest  are  governed  by  non-linear  operators.  However,  a  linear  analysis  (e.g.,  stability 
analysis  of  explicit  methods,  modified  equation,  convergence  rates  for  iterative  methods, 
etc.)  is  often  used  successfully  to  make  the  problem  tractable  and  to  provide  insight  into 
the  problem  of  interest.  Equation  (A.3)  shows  that  the  modeling  error  8Sm  is  governed  by 
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the  same  operator  as  the  solution  M,  but  with  a  source  term  and  differences  in  the  IC  and 
BC  due  to  errors  in  modeling  the  true  operator  and  IC  and  BC  functions. 

The  analytical  approach  solves  equation  (A. 2)  exactly  and  is  thus  limited  to  simple 
fluid  mechanics  problems.  The  simulation  approach  solves  equation  (A.2)  approximately 
using  numerical  methods  and  thus  introduces  numerical  errors.  The  continuous  IBVP  is 
reduced  to  a  discrete  IBVP,  i.e.,  algebraic  equations  using  spatial  and  temporal 
discretization  techniques  such  as  finite  difference,  volume,  and  element  methods  resulting 
in  numerical  errors  due  to  spatial  Ax,  temporal  At,  and  other  step  sizes  Ax.  (i.e.,  the 

numerical  error  is  zero  when  the  step  sizes  are  zero).  The  discrete  IBVP  is  defined  by  a 
discrete  operator  with  discrete  IC  and  BC 

A(S)  =  r7 

IC:S(x,t  =  O)  =  0N(x)  (A.  4) 

BC:S(xb,  t)  =  '»N(t) 

where  the  source  term  T/  is  the  residual  imbalance  of  the  algebraic  equations  due  to  the 
use  of  implicit  methods.  If  explicit  methods  are  used,  iterative  errors  do  not  exist  and  T/  = 
0.  Equation  (A.4)  is  solved  on  a  computer  through  a  set  of  programming  instructions  (i.e., 
a  CFD  computer  code)  to  provide  the  simulation  prediction  S.  Program  execution 
requires  specification  of  various  input  parameters,  including  step  size  distributions. 

Numerical  errors  can  be  defined  and  evaluated  by  transforming  the  discrete  IBVP 
back  to  a  continuous  IBVP.  This  is  accomplished  by  representing  A  as  a  generalized 
Taylor  series  about  a  numerical  benchmark  Sc  (solution  with  zero  step  sizes)  in  terms  of 
step  sizes  Ax . 


S  =  SC 


-it 


(Ax,.)'  &S 

/!  dxj 


(A. 5) 


where  j=l,  J  is  used  to  represent  various  step  sizes  introduced  in  discretization  of  the 
continuous  PDE,  IC,  and  BC  (spatial  AxG  =Ax,  temporal  Axf  =  At ,  and  other  Axj). 
Substituting  expansion  (A.5)  into  equation  (A.4)  and  rearranging  gives  the  modified 
equation  that  is  actually  solved  when  discretization  techniques  are  applied  to  an  IBVP 
(Anderson  et  al.,  1984) 

LModiflediS)  =  LM(S)  =  TN 

IC:S(x,t  =  0)  =  GModifwd(x)  (A. 6) 

BC:S(xB,t)  =  HModifwd(t) 

where  the  source  term  is  given  by 

r„=r,+£r,  (a.7) 

The  summation  term  in  equation  (A.7)  represents  the  truncation  errors  due  to  differences 
between  the  continuous  and  discrete  PDE.  Spatial  and  temporal  truncation  error  terms  for 
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typical  finite  difference  and  volume  methods  are  in  the  form  of  a  power  series  expansion  in 
step  sizes  A Xj 


rj =E(Av)‘ 


a 


(0 


i= 1 


(A. 8) 


where  the  superscript  (/)  is  used  to  indicate  variables  in  the  i,h  term  of  the  expansion, 


,dlS 


ay  =  / ( — -)  contains  solution  derivatives  with  respect  the  xj  and  are  independent  of  the 

dxj 

step  sizes  Ax.,  and  pl‘ '  is  the  rate  of  reduction  of  the  truncation  error  terms  with 
refinement  of  Ax .  (i.e.,  order  of  accuracy).  The  modified  equation  (A. 6)  recovers  the 


modeled  operator  LM\  however,  it  operates  on  the  simulation  prediction  S  instead  of  the 
exact  solution  to  the  modeled  equations  M.  Thus,  the  source  term  rv  causes  the 
simulation  prediction  S  to  differ  from  the  exact  solution  to  the  modeled  equation  M. 


Subtracting  equations  (A. 2)  and  (A. 6)  gives  the  IB  VP  that  governs  the  simulation 
numerical  error  Sw  =  S  -  M  (i.e.,  the  difference  between  the  simulation  and  modeled 
values)  (Ferziger  1993;  Roache  1998) 


lm  (S  -  M)  =  L„  (5  „ )  =  r„  =  r,  +  £  r, 

M 

IC :  8  „  (x,  t  =  0)  =  (x)  -  G„  (or)  (  A.  9) 

BC  :  8  SN  ( XB  ,  t)  =  H  Madifwd  (0  —  H M  (t) 


Thus,  the  iterative  and  truncation  error  terms  also  act  as  source  terms  for  numerical 
errors  8 sn  in  the  solution  S.  If  there  are  no  iterative  errors,  the  source  term  V N  (and 
thus  8 sn)  is  zero  when  either  the  truncation  error  is  zero  (e.g.,  spectral  methods)  or  step 
size  Ax  .  is  zero.  If  the  exact  form  of  the  truncation  error  terms  [equation  (A. 8)]  for  a 

discretization  technique  is  known,  equation  (A.  9)  can  be  solved  numerically  to  give  S.sv. 
Such  methods  can  be  classified  as  single  step-size  error  estimation  methods  (e.g., 
Shimazaki  et  al.,  1993). 

Rewriting  the  PDE  for  the  numerical  error  in  equation  (A. 9)  with  the  source  term 
r N  expanded  gives 

i„(5sv)  =  r,  +EE<A x/'"af  (A.  10) 

M  «'= 1 

Since  the  term  (Av;  )p'  is  independent  of  the  continuous  operator  LM,  the  numerical 

error  Saw  [i.e.,  the  solution  to  equation  (A.  10)]  reduces  at  the  same  rate  as  the  source  term 
Twin  equation  (A.  10)  so  that  the  solution  is  of  the  form 

SW=S/+ESy  (A-U) 

/=i 

where 
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(A.  12) 


8.  -I>V)  '  gj 


(0 


1= 1 


0/ 

and  g(,)  =/( — -)  is  the  “grid”  function  which  contains  continuous  solution  derivatives. 

chj 

The  form  of  equation  (A.  11)  can  be  verified  by  substitution  of  equation  (A.ll)  with 
(A.  12)  into  (A.  10),  which  gives 

£  „  [8  ,+tt  < Al  /  y'“  g?  ]  =  LM  (8, )  +  E  £  K  [(A* , ) gf  ]  (A.  1 3) 

/= 1  i= 1  7=1  i=l 

If  step  size  Ax .  and  order  of  accuracy/?/  are  constant  and  independent  of  LM,  the  last  term 
in  equation  (A.  13)  can  be  rewritten 

Lu  [8 ,  + 1,  E  (A* ,  )"rgf  1  =  Lu  (5 , )  +  £  £  K '  L><  bfj  (A.I4) 

y=l  z=l  7=1  i=l 

Comparison  with  equation  (A.  10)  gives  the  following  equations 

LM@i)  =  ri  (A.  15) 

^[gf]  =  af  (A.  16) 


The  solution  to  equation  (A.  16)  for  g{‘ 1  is  independent  of  step  size  A x  .  and  order  of 
accuracy/?/.  Thus,  the  form  of  the  numerical  error  bSN  assumed  in  equation  (A.ll)  involve 
products  of  (Ax-,  )p'  and  functions  that  are  independent  of  step  size  and  order  of 

accuracy.  This  verifies  that  the  numerical  error  S.sV  reduces  as  (Ax/  )Pj  . 


Finally,  the  IBVP  that  governs  the  simulation  error  8V  is  obtained  by  adding 
equations  (A. 3)  and  (A. 9) 

LM(S-T)  =  LM$S)  =  rN+rM 

IC  :  5  s  (x,0)  =  GModified  (x)  -  Gt  (x)  (A.  17) 

BC  :  8  s  ( xB ,  t )  =  H Modified  it)  —  Ht  ( t ) 

where  simulation  error  is  defined  as 

ds=S-T  =  8SN+bSM  (A.  18) 

Equation  (A.  18)  provides  the  desired  expression  for  the  simulation  error  in  terms  of  the 
simulation  modeling  and  numerical  errors.  It  shows  that  the  simulation  modeling  and 
numerical  errors  are  additive  subject  to  the  assumption  that  LM  is  a  linear  operator. 
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Appendix  B.  Generalized  Richardson  Extrapolation 

For  the  converging  condition  (i)  in  equation  (16),  generalized  RE  is  used  to  estimate 
Uk  or  8  *  and  Uk  .  RE  is  generalized  for  J  input  parameters  and  for  use  of  correction 

factors  to  account  for  the  effects  of  higher-order  terms  and  defining  and  estimating  errors 
and  uncertainties,  as  summarized  in  Section  3.2.3.  This  appendix  provides  a  detailed 
description. 

Generalized  RE  begins  with  equation  (14).  The  error  terms  on  the  right-hand-side  of 
equation  (14)  are  of  known  form  (i.e.,  power  series  expansion  in  Axk )  based  on  analysis 
of  the  modified  (A. 6)  and  numerical  error  (A. 9)  equations,  as  shown  in  Appendix  A 
equation  (A.  12),  which  is  written  below  as  a  finite  sum  (i.e.,  error  estimate)  and  for  the  kth 
parameter  and  mth  solution 

=£(^i.)'iAr  (B.i) 

i= 1 

n  =  number  of  terms  retained  in  the  power  series,  powers  /^/'correspond  to  order  of 

accuracy  (for  the  ith  term),  and  g)‘ 1  arc  referred  to  as  “grid”  functions  which  are  a 

function  of  various  orders  and  combinations  of  derivatives  of  S  with  respect  to  jc*. 
Substituting  equation  (20)  into  equation  (14)  results  in 

k.=Sc+tl(tel,y°g>‘,+  XX  (B.2) 

(=1  j=hj*k 

Subtraction  of  multiple  solutions  where  input  parameter  Axk  is  uniformly  refined 
eliminates  the  8*  terms  in  equation  (B.2)  since  8*  is  independent  of  Ax',  and  provides 

equations  for  Sc,  p)n ,  and  gk] .  This  assumes  pk  and  gyk  are  also  independent  of  Ax',  . 

Since  each  term  (i)  contains  2  unknowns,  m=2n+l  solutions  are  required  to  estimate  the 
numerical  benchmark  Sc  and  the  first  n  terms  in  the  expansion  in  equation  (B.2)  (i.e.,  for 
n=  1,  m= 3  and  for  n= 2,  m= 5,  etc).  The  accuracy  of  the  estimates  depends  on  how  many 
terms  are  retained  in  equation  (B.I),  the  magnitude  (importance)  of  the  higher-order 
terms,  and  the  validity  of  the  assumption  that  p[n  and  g(k  are  independent  of  Axk  .  For 
sufficiently  small  Axk ,  the  solutions  are  in  the  asymptotic  range  such  that  higher-order 
terms  are  negligible  and  the  assumption  that  pyk  and  gk  are  independent  of  Axk  is  valid. 

However,  achieving  the  asymptotic  range  for  practical  geometry  and  conditions  is  usually 
not  possible  and  m> 3  is  undesirable  from  a  resources  point  of  view;  therefore,  methods  are 
needed  to  account  for  effects  of  higher-order  terms  for  practical  application  of  RE. 
Additionally,  methods  may  be  needed  to  account  for  possible  dependence  of  pk]  and  g*/' 

on  Ax',  ,  although  not  addressed  herein.  Usually  8  k  is  estimated  for  the  finest  value  of  the 
input  parameter,  i.e.,  8  *=8**  corresponding  to  the  finest  solution  Sk  .  RE  can  be 
classified  as  a  multiple  step-size  error  estimation  method. 
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If  only  the  leading  term  (n= 1)  in  equation  (B.l)  is  estimated,  three  solutions  are 
required  and  can  be  written  from  equation  (B.2) 

sk{  =  Sc  +  (ax,  Y"  gf  +  £8-  (B.3) 

j=U*k 

(  \p(l)  J 

K  =Sc+(rlAx,y  g">  +  £5’  (B.4) 

■S,  =  Sc  +  (x>,  k' gf  +  £8  *  (B.5) 

j=l,i*k 


Equations  (B.3)-(B.5)  provide  three  equations  for  the  three  unknowns  (Sc,  p[' ' ,  and 
g{n ).  The  order  of  accuracy  and  “grid”  function  are  found  by  computing  the  solution 
changes  e21  =  Sk_  —Sk  and  E32i  =  Sk  -  St  from  equations  (B.3)  -  (B.5) 


=4-Y  =kk  «;■><<■/"  -i) 

(B.6) 

=rf  Kb"  ('i'”  “I) 

(B.7) 

The  estimate  of  the  error  8^ 


~  8  ' '  ’  =  (Ax^  )Pk  gf*  is  obtained  from  equation  (B.6) 


8 


HD  _ 

RE, 


'2U 


n(i) 

fr/'  -1) 


(B.8) 


where  8 ^ 1 '  is  an  estimate  of  the  first  term  of  the  expansion  in  equation  (B.l)  using  RE. 

The  order  of  accuracy  p)u  is  obtained  by  eliminating  the  term  (Ax^  )Pt  gt)(rkk  “!)  from 
equation  (B.6)  and  (B.7) 


Pk]  = 


fr(^32t  ^2lk  ) 
ln(/-J 


(B.9) 


If  order  of  accuracy  is  assumed  known  (e.g.,  from  the  modified  equation  or  from  grid 
refinement  tests  for  simple  geometry  using  similar  grid  expansion)  only  two  solutions  are 
required  to  obtain  an  estimate  of  the  leading  term  in  the  power  series  expansion  equation 
(B.l).  However,  a  minimum  of  three  solutions  is  required  to  establish  convergence  with 
refinement  of  input  parameter. 

An  estimate  using  the  first  two  terms  (n= 2)  in  equation  (B.l)  can  be  obtained  from 
five  solutions 


Sc  +(Axt  ) 


gf  + 


(a**,  ) 


pI2' 


gk] 


+  is; 

j=U*k 


(B.10) 
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Sc 

+ 

(rA v 

\p? 

g? 

+  ( 

Pk^k,] 

|  Pi2) 

g?  ' 

+ 

tK 

(B.l  1) 

j 

r=l  ,i*k 

/ 

\n(2) 

J 

Sc 

+ 

V rfakl 

Y‘ 

g? 

+ 

l rk^kx 

Y‘ 

g f 

+ 

Is;, 

(B.12) 

j=\,i±k 

/ 

/ 

\n(  2) 

J 

Sc 

+ 

Vk  A**, 

)"* 

g? 

+ 

\P  A**, 

)'* 

g? 

+ 

Is;, 

(B.13) 

j=U*k 

/ 

/ 

J 

Sc 

+ 

VV  A**. 

)'* 

gf 

+ 

Yk^k, 

Y‘ 

gf 

+ 

Is;, 

(B.14) 

j=U*k 

The  orders  of  accuracy  and  the  grid  functions  are  obtained  by  computing  the  four  solution 
changes  e2l(  =  Ski  -  Sk< ,  £32;  =  Sk  -S^,  £43t  =  Skf  -Sk},  and  £54t  =5*.  -Ski  which  gives 
four  equations  for  the  four  unknowns,  pk] ,  p',2' ,  g)''  and  g[2) .  Upon  solution,  the  four 
unknowns  are  used  to  give  an  estimate  of  the  first  and  second  terms  in  equation  (B.l) 


„Pk 


*(2) 


£2,  £ 


*Pk 


32, 


£  21 ,  e 


32, 


’RE, 


Pk 


■1) 


Pk 


(B.15) 


'  k  ’  k  A' A-  V  Vk  ’k  >\’k  "1) 

where  8  '[2)  is  an  estimate  of  the  first  two  terms  of  the  expansion  in  equation  (B.l)  using 

RE.  The  orders  of  accuracy  of  the  first  and  second  terms  in  the  expansion  pk]  and  p{2) 
are  given  by 


n(„_ln[a, -pj 

Pk 


In (r.) 


Pk 


(2) 


’k> 

InK  +  Pj 


(B.16) 


ln(ft) 


where 


ak  =ak/bk 


ak  ~  e21i£54i 


'E32te43t 


4 

■  3e  2 


E21tE434  E32t  , 


32tE43t  +4£21t£43j,  +  ^E32tE54t 


8e21te32te43te54t 


+  £21  E54t 


Two  conditions  are  required  to  obtain  estimates  of  the  orders  of  accuracy  from  equation 
(B.16):  (i)  ak  ±(3t  >1;  and  (ii)  ck  >  0 .  Condition  (i)  is  be  satisfied  if  the  solutions  are 

monotonically  convergence  while  condition  (ii)  is  satisfied  if  the  solutions  are  sufficiently 
close  to  the  asymptotic  range.  If  the  orders  of  accuracy  are  assumed  known,  only  three 
solutions  are  required  to  estimate  the  first  two  terms  in  the  power  series  expansion  using 
equation  (B.15). 
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Appendix  C.  Analytical  Benchmarks 

Analytical  benchmarks  can  be  defined  as  the  truth  and  are  useful  in  development 
and  confirmation  of  verification  procedures  and  methodology  and  in  code  development, 
but  can  not  be  used  for  validation  and  are  restricted  to  simple  equations.  Results  were 
obtained  for  two  analytical  benchmarks:  one-dimensional  (ID)  wave  and  two-dimensional 
(2D)  Laplace  equations.  The  results  for  the  2D  Laplace  equation  were  qualitatively  similar 
to  those  for  the  ID  wave  equation,  which  are  presented  in  this  appendix. 

The  IBVP  and  solutions  for  the  true,  model,  and  analytical  benchmarks  are 
equivalent  such  that  the  modeling  error  is  zero  and  the  only  simulation  error  is  the 
simulation  numerical  error.  Simulation  results  (for  the  ID  wave  and  2D  Laplace 
equations)  are  compared  to  analytical  benchmark  solutions  to  determine  the  exact 
simulation  numerical  error  and  evaluate  both  single  and  multiple  step-size  error  estimation 
methods.  In  the  latter  case,  generalized  RE  is  used  and  the  role  of  higher-order  terms  in 
the  power  series  expansion  of  the  simulation  numerical  error  [equation  (B.l)]  is  assessed 
by  comparing  estimates  of  the  leading  term  using  three  grids  to  those  for  the  first  two 
terms  using  five  grids.  Correction  factors  are  derived  to  account  for  the  effects  of  the 
higher-order  terms  and  to  define  the  uncertainty  in  the  error  estimate. 


C.l  Verification  Using  Analytical  Benchmarks 

By  definition,  the  IBVP  for  the  true,  model,  and  analytical  benchmark  solutions  are 
equivalent 

Lt(T)=Lm(M)  =  La(A)  =  0 

1C  :  T(x,  t  =  0)  =  M(x,t  =  0)  =  A(x, t  =  0)  =  Gr ( x )  =  Gu  (x)  =  GA (x)  (C.l) 

BC  :T(xB,t)  =  M(xB,t)  =  A(xB,t)  =  HT(t )  =  HM  ( t )  =  H A(t) 


Therefore, 

T  =  M  =  A 

(C.2) 

and 

s»=o 

(C.3) 

The  simulation  error  and  uncertainty  are  given  by 

Ss=S-A  =  Ss„ 

(C.4) 

ul=ul 

(C.5) 

and  the  corrected  simulation  error  and  uncertainty  are  given  by 

85c  =  $C  ~  A  —  &SN 

(C.6) 

u2  =u 2 

U  Sc  U  SCN 

(C.7) 

Simulations  are  verified  if 
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(C.8) 


\E\  =  \A  -  S\  <  USN 

and  corrected  simulations  are  verified  if 

\Ec\  =  \A-Sc\<UScN  (C.9) 


C.2  IBVP  and  Analytical  and  Numerical  Solutions  for  ID  Wave  Equation 

The  ID  wave  equation  is  called  a  model  equation,  as  it  models  the  behavior  of 
more  complicated  (nonlinear)  PDE.  A  simplified  form  is  the  first-order  linear  convection 
equation  with  IBVP 


c)A  dA 

- I-C - 

dt  dx 


L,M)  =  —  +  c—  =  0 


1C :  A(x,0)=  A0  exp 


M 

B 


BC :  A(-°°,t)  =  0 


(C.10) 


The  initial  condition  is  prescribed  by  a  Gaussian  function  centered  at  x  =  0.0  with  A0  =  1 
and  B  =  0.005,  the  boundary  condition  far  upstream  is  zero,  and  c  is  the  wave  speed, 
which  is  set  to  unity.  The  computational  domain  is  defined  as  -1  <  x  <  2  . 

The  exact  solution  to  equation  (C.10)  is 


A(x,t)  =  A)  exP 


(C.ll) 


Figure  C.l  shows  the  initial  condition  and  the  exact  solution  at  t  =  1. 

Two  discretization  techniques  are  studied:  (i)  a  first-order  (Euler)  explicit  method; 
and  (ii)  a  second-order  implicit  method. 

For  the  Euler  explicit  method,  the  discrete  operator  in  equation  (A.4)  is  given  by 


m+\ 


A(S)  =  - 


Att 


- s r  s"  -  s" ,  . 

-  +  c— - —  =  0 


A Xr 


(C.12) 


where  n+1  and  n  denote  the  new  and  current  time  levels,  respectively.  The  modified 
equation  (A. 6)  and  simulation  numerical  error  equation  (A. 9)  are  given  by 


ds  ds  _ 

*+cah=r« 

(C.13) 

,  dbSN  _  T- 
dt  dx  N 

(C.14) 

with  source  terms  [equation  (A. 8)] 

f  rS  \ 
r  =  At  _ 

V  1  ) 


+  At7 


+  0[(Atg  )2 ,  (Avr  )2  ] 


(C.15) 
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For  equations  (C.17)  and  (C.18),  the  IC  and  BC  are  the  same  as  for  equations  (C.13)  and 
(C.14),  respectively. 

The  form  of  the  solution  to  equation  (C.14)  and  (C.17)  was  given  in  Appendix  A 
by  equation  (A.  11)  and  (A.  12),  as  a  power  series  expansion  in  A xj.  For  the  Euler  explicit 
method,  the  form  is  given  by 

85jv  =  5g  +5r  =AxgSg’  +A xTgf  +  0[{Axg  )2 ,  (Axr  )2  ]  (C.20) 

For  the  second-order  implicit  method,  the  form  is  given  by 

5^  =SG  +5r  =(AxG)2g®  +(A xTfgf  +  0[(Axg  )4 ,  (Axr  )4  ]  (C.21) 

Results  were  obtained  for  the  numerical  solution  of  equation  (C.12)  and  C.16) 
using  ten  grids  (Table  C.l)  and  with  two  values  of  CFL  =  cAt  /  Av  =  (0. 1,0.5) .  The 
solutions  were  monotonically  convergent  for  all  ten  grids  and  both  CFL  based  on  the 
convergence  ratio  R  [equation  (15)]  defined  with  the  ratio  of  the  L2  norm  of  solution 

changes  e32w  and  e21w  .  Figure  C.2  compares  simulation  S  to  analytical  benchmark  A 

for  t=  1  (along  with  single  and  multiple  step  size  error  estimates  to  be  discussed  later).  The 
inherent  deficiencies  of  the  two  methods  are  apparent.  The  first-order  method  displays 
dissipation  errors  due  to  even  simulation  derivatives  in  equation  (C.15),  which  reduce  at  a 
first-order  rate,  whereas  the  second-order  implicit  method  displays  dispersion  errors  due 
to  odd  simulation  derivatives  in  equation  (C.19),  which  reduce  at  a  second-order  rate. 


C.3  Single  Step  Size  Error  Estimation  Method 

Single  step  size  error  estimation  methods  are  based  on  solution  of  the  IB  VP  for  the 
simulation  numerical  error  &Sn,  as  given  by  equation  (A.9).  This  provides  an  estimate  of 
the  numerical  error  for  a  single  step  size.  Two  step  sizes  are  required  to  evaluate 
convergence  with  respect  to  input  parameter.  Such  methods  require  fewer  solutions  than 
multiple  step-size  error  estimation  methods;  however,  there  are  several  obstacles  for 
practical  problems.  Derivation  of  the  modified  equation  (A. 6)  is  necessary  in  order  to 
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define  the  truncation  error  terms  [equation  (A. 8)],  which  are  the  source  terms  for  equation 
(A.9).  This  may  be  difficult  or  not  possible  depending  on  the  complexity  and  type  of 
discretization  technique  used.  The  coefficients  in  the  source  term  are  functions  of  higher- 
order  solution  derivatives.  Higher-order  discretization  techniques  with  associated 
increased  numerical  instabilities  must  be  used  to  discretize  the  numerical  error  equation 
(A.9)  than  those  in  used  in  the  original  IB  VP  equation  (A.4)  such  that  the  truncation  error 
terms  for  the  discrete  form  of  equation  (A.9)  are  higher  order  than  those  for  equation 
(A.4).  Also,  additional  programming,  memory,  and  computer  time  are  required  to  include 
solution  of  the  simulation  numerical  error  equation. 

Results  were  obtained  for  the  numerical  solution  of  the  numerical  error  equation 
(C.14)  and  (C.18)  for  grids  6-10  of  Table  C.l  and  with  CFL= 0.1.  Fourth-order  spatial  and 
third-order  time  discretization  techniques  were  used.  Figure  C.2  compares  the  exact 
comparison  error  E=A-S  to  the  single  step  size  error  estimate  for  8,w.  The  results  show 
that  single  step  size  error  estimates  are  accurate  even  for  the  first-order  method  and 
coarsest  grid. 


C.4  Multiple  Step  Size  Error  Estimation  Method 

Multiple  step  size  error  estimation  methods  are  based  on  generalized  RE,  as 
described  in  Appendix  B.  The  total  true  numerical  error  (i.e.,  grid  size  and  time  step)  can 
be  computed  since  the  exact  solution  is  known  for  the  analytical  benchmark.  However,  the 
exact  grid  size  or  time  step  error  cannot  be  computed  separately.  As  such,  a  combined 
grid  and  time  step  study  was  conducted  with  CFL= 0.5  for  all  ten  grids  to  directly  compare 
the  true  error  to  estimates  from  RE.  For  the  combined  grid  and  time  step  study,  overall 
order  of  accuracy  (i.e.,  spatial  and  temporal)  is  estimated  and  the  subscript  SN  is  used  to 
denote  an  estimate  of  total  simulation  error. 


The  role  of  higher-order  terms  in  the  power  series  expansion  of  the  simulation 
numerical  error  is  assessed  by  comparing  estimates  of  the  leading  term  using  three  grids  to 
those  for  the  first  two  terms  using  five  grids.  To  avoid  problems  associated  with 
pointwise  calculation  of  order  of  accuracy  discussed  in  Section  3.2.3,  the  orders  of 
accuracy  are  defined  using  the  L2  norm  of  the  solution  changes.  Order  of  accuracy  of 
the  first  term  in  the  error  expansion  is  given  by 
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which  is  used  to  provide  a  pointwise  error  estimate  8  ^  of  the  leading  term  in  the  error 
expansion 
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Similarly,  the  equations  for  the  first  two  terms  in  the  error  expansion  are 


(C.23) 
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Recall  from  Appendix  B  that  the  conditions  for  applying  equation  (C.24)  were  that 
(i)  aSN  ±  Psw  >  1  (related  to  monotonic  convergence)  and  (ii)  cSN>  0  (related  to  the 
solutions  being  in  the  asymptotic  range).  It  was  found  that  condition  (ii)  was  not  satisfied 
for  grid  sizes  Ax  >  3 . 1 24x  1  0~4  for  the  first-order  scheme. 

Tables  C.2  and  C.3  show  the  orders  of  accuracy  for  the  first-order  and  second- 
order  methods,  respectively.  For  the  first-order  method,  the  three-grid  estimate 
approaches  the  theoretical  rate  p\l '  =  1  from  below  as  the  grid  is  refined,  whereas  the  five- 
grid  estimate  approaches  the  theoretical  rates  p)'h '  =  1  and  p]l'  =2  from  above  and  below, 
respectively,  as  the  grid  is  refined.  For  the  second-order  method,  the  solutions  are  in  the 
asymptotic  range  even  on  the  coarsest  grids,  although  p[2)  is  larger  than  p)l'  =4. 


Figure  C.2  compares  the  exact  comparison  error  E=A-S  to  the  three-grid  error 
estimate  and  the  single  step  size  error  estimate  for  bSN-  For  the  first-order  method,  the 
three-grid  estimate  is  relatively  poor  especially  for  the  coarser  grids,  whereas  for  the 
second-order  method  the  three-grid  estimate  is  close  to  both  E  and  6Sn-  Figure  C.3 
compares  E  to  both  the  three-  and  five-grid  error  estimates.  For  the  first-order  method, 
the  three-grid  estimate  is  less  accurate  than  the  five-grid  estimate  especially  for  the  coarser 
grids,  whereas  for  the  second-order  method  both  the  three  and  five  grid  estimates  are 
accurate.  The  results  show  that  the  higher-order  terms  are  more  important  for  lower- 
order  methods  on  coarser  grids. 
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C.4  Correction  Factors,  Uncertainties,  and  Verification 

Methods  are  needed  to  account  for  the  effects  of  higher-order  terms  for  practical 
application  of  RE  when  solutions  are  outside  the  asymptotic  range.  Figure  C.4a  and  C.5a 
compare  the  true  error  E  to  the  three-grid  error  estimate  8^n  vs.  step  size  at  one  spatial 

location  (x=  1  for  the  first-order  method  and  x=l.l  for  the  second-order  scheme  since 
maximums  of  numerical  error  occur  there)  and  for  the  first-  and  second-order  methods, 
respectively.  The  three-grid  estimate  accurately  estimates  the  true  error  E  for  smaller  step 
sizes,  but  over  predicts  E  for  larger  step  sizes.  Closer  examination  reveals  that  the  reason 
equation  (C.23)  over  estimates  the  error  is  due  to  the  fact  that  equation  (C.22)  under 
estimates  the  order  of  accuracy,  as  also  shown  in  the  figures  and  previously  in  Table  C.2. 
Therefore,  one  approach  is  to  correct  the  three-grid  estimate  by  a  multiplication  correction 
factor,  which  accounts  for  this  deficiency,  i.e. 

K=CK  (C26) 

Two  definitions  for  C*  were  investigated.  The  first  is  based  on  equation  (B.8)  substituted 
for  the  left-hand  side  of  equation  (C.26)  and  solving  for  Ca,  but  replacing  p[u  from 
equation  (B.9)  by  an  improved  estimate  pk  based  on  the  modified  equation  pk  (or  in  the 
general  case  on  solutions  for  simplified  geometry  with  similar  step  size  expansion) 
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Similarly,  the  second  is  based  on  equation  (B.15),  but  replacing  p(k  and  p[2)  from 
equation  (B.16)  by  improved  estimates  pk  and  qk 
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Figure  C.4a  and  C.5a  also  compare  E  with  error  estimates  based  on  equation  (C.26)  with 
equation  (C.27)  and  (C.28).  Both  estimates  are  closer  to  E  than  the  uncorrected  three 
grid  estimate  8  *|l) ,  but  for  coarser  grids  C[l)  is  somewhat  too  small  and  C)2'  is  slightly 
too  large.  Figure  C.4b  and  C.5b  show  the  same  trends,  but  directly  compare  the  exact 
correction  factor E /S*^  to  equation  (C.27)  and  (C.28).  In  this  case,  Ca<1  indicates  that 
the  leading-order  term  over  predicts  (higher-order  terms  net  negative)  the  error. 
However,  for  the  general  case,  Ca  is  equally  likely  to  be  <  1  or  >  1  depending  whether  the 
order  of  accuracy  is  approached  from  below  or  above,  respectively.  Ck  >1  indicates  that 
the  leading-order  term  under  predicts  (higher-order  terms  net  positive)  the  error.  Thus, 
for  the  general  case  the  correction  to  the  leading-term  error  estimate  is  equally  likely  to  be 
positive  or  negative  and  can  be  used  to  define  the  simulation  numerical  uncertainty. 

Equation  (C.26)  is  used  to  estimate  Uk  or  8  *  and  Uk  depending  on  how  close  the 
solutions  are  to  the  asymptotic  range  (i.e.,  how  close  Ck  is  to  1)  and  one’s  confidence  in 
equation  (C.26).  There  are  many  reasons  for  lack  of  confidence,  especially  for  complex 
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three-dimensional  flows.  Often,  pointwise  results  are  not  uniformly  convergent  over  all 
grid  points  (i.e.,  locally  oscillatory  or  even  divergent). 

For  Ck  sufficiently  less  than  or  greater  than  1  and  lacking  confidence,  Uk  is 

estimated,  but  not  8  *  and  Uk  .  Figure  C.4d  and  C.5d  show  that  equation  (C.26)  can  be 

used  to  estimate  the  uncertainty  by  bounding  the  error  by  the  sum  of  the  absolute  value  of 
the  corrected  estimate  from  RE  and  the  absolute  value  of  the  amount  of  the  correction 
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For  Ck  sufficiently  close  to  1  and  having  confidence,  8  *  and  Uk  are  estimated. 
Equation  (C.26)  is  used  to  estimate  the  error  SJ,  which  can  then  also  be  used  in  the 
calculation  of  Sc  [in  equation  (10)].  Figure  C.4c  and  C.5c  show  that  uncertainty  in  the 
error  estimate  can  be  based  on  the  amount  of  the  correction 
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Note  that  in  the  limit  of  the  asymptotic  range,  Ck  =1,  8k  =  8k  =  8ffi  ,  and  Uk  =0. 

1  h  c 

Uncertainty  estimates  enable  a  quantitative  measure  of  verification  for  analytical 
benchmarks.  A  simulation  is  verified  if  equation  (C.8)  is  satisfied  and  corrected 
simulations  are  verified  if  (C.9)  is  satisfied.  Figure  C.4c,d  and  C.5c,d  indicate  that  the 
present  solutions  are  verified  at  the  chosen  spatial  location  and  at  the  levels  ( Uk , 
Uk  )=(15%,  7.5%)  and  (1.2%,  .2%)  for  the  first-  and  second-order  methods,  respectively. 
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Grid 


number  points 


Ax 


1 

38401 

7.8125  x  10"5 

2 

19201 

1.5625  x  10"4 

3 

9601 

3.125  x  10"4 

4 

4801 

6.25  x  10-3 

5 

2401 

1.25  x  10"3 

6 

1201 

2.5  x  10-3 

7 

601 

5  x  10-3 

8 

301 

1  x  10-2 

9 

151 

2  x  10-2 

10 

76 

4  x  10-2 

Table  C.l.  Grids 


Ax 

Estimate 

P? 

P? 

6.25xl0"4 

3  grid 

0.94 

- 

5  grid 

1.06 

1.46 

3.125xl0"4 

3  grid 

0.97 

- 

5  grid 

1.01 

1.73 

Table  C.2.  Orders  of  accuracy  for  first-order  method 


Ax 

Estimate 

P? 

P? 

2.5xl0"3 

3  grid 

2.00 

- 

5  grid 

2.00 

4.54 

Table  C.3.  Orders  of  accuracy  for  second-order  method 
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Error  _  Solution, 


Figure  C.l.  Initial  condition  and  exact  solution  for  the  ID  Wave  Equation 


Figure  C.2.  Numerical  solution  S,  true  error  E,  single  grid  estimate  6  *v  ,  and  three-grid 

estimate  from  RE  for  ID  wave  equation  at  t=  1:  (a),  (b)  Euler  explicit  scheme 
(first-order)  and  (c),  (d)  second-order  implicit  scheme. 
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Error 


Figure  C.3.  Comparison  of  the  true  numerical  error  E  to  error  estimates  from  RE  using 
three  and  five  grids  for  ID  wave  equation  at  t=  1 : 

(a)  Euler  explicit  scheme  (first-order)  and  (b)  second-order  implicit  scheme. 


Figure  C.4.  Verification  results  for  first-order  numerical  solution  of  ID  wave  equation,  (a) 
comparison  of  true  error  A-S  to  estimates  from  RE,  (b)  correction  factor,  and 
(c)  comparison  of  ASC  and  ±UScn,  and  (d)  comparison  of  A-S  and  ±USn- 
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-0.015 


Figure  C.5.  Verification  results  for  second-order  numerical  solution  of  ID  wave  equation, 
(a)  comparison  of  true  error  A-S  to  estimates  from  RE,  (b)  correction  factor, 
and  (c)  comparison  of  A-Sc  and  ±USCn,  and  (d)  comparison  of  A-S  and  ±USn- 
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